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Introduction 

by 

Larry L. Schumaker 

On may 17-19, 1982, a workshop entitled "Surface Fitting" was held 
at Texas A&M University, College Station. The puroose of the workshop 
was to bring together leading experts from academia, industry, and 
government laboratories for an exchange of views and a discussion of the 
"state of the art." For a list of participants see pages 5-6. 

The workshop began with an overview by R. P. Heydorn of NASA/Johnson 
Space Center, Houston, Texas. The purpose of the overview was to acquaint 
the participants with some mathematical/statistical problems within the 
AgRISTARS Program which may be amenable to investigations involving the 
use of surface-fitting techniques. In order to establish a framework 
for the workshop, Larry Schumaker presented a general survey of surface 
fitting and contouring in which he touched on a variety of local and 
global methods for both interpolation and approximation. 

The program for the workshop included six invited lecturers (see the 
program on pages 3-5). Charles Lawson discussed the construction of a 
triangular grid on the sphere, and the computation of corresponding C^ 
surfaces. R. E. Barnhill dealt with several schemes based on patches and 
blending. Rolland Hardy lectured on the multiquadric surfaces which he 
invented. Thomas Foley considered a three stage procedure which proceeds 
from scattered data to grid values using local least squares, then to a 
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bicubic B-spline interpolant, and finally uses Shepard's method to 
obtain an interpolant of the original data. Douglas Bates discussed 
smoothing splines and the method of generalized cross validation, with 
particular emphasis on computational methods. Rosemary Chang's talk 
dealt with several practical problems arising in Engineering. 

In addition to the formal lectures, the program included a panel 
discussion involving everyone. We believe that the workshop gave all 
participants-- the theoreticians, the practitioners, and the consumers— 
a better understanding of what methods and software are available, and of 
what needs to be done in the future. 

Written versions of the lectures are included in this document. 

The talks elicited a great deal of discussion which we have not attempted 
to reproduce here. Finally, this document includes a computerized 
bibliography of surface fitting papers which Larry Schumaker has 
assembled at Texas A&M University. Additions and corrections to this 
bibliography would be greatly appreciated. 
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Coffee break 

- 

11:00 - 12:00 
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Dinner at Larry Guseman's cabin 


Wednesday, May 19 

8:15 - 8:30 Coffee & doughnuts 
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to NASA 


10:00 - 10:30 Coffee Break 
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FITTING SURFACES TO SCATTERED DATA 
Larry L. Schumaker 

This paper is a survey of a variety of numerical methods 
for fitting a function to data given at a set of points scat- 
tered throughout a domain in the plane. We discuss four 
classes of methods: (1) global interpolation, (2) local inter- 

polation, (3) global approximation, and (4) local approximation. 
We also discuss two-stage methods and contouring. The surfaces 
constructed will include polynomials, spline functions, and ra- 
tional functions, among others. 

1. Introduction 

Out aim is to survey methods for solving the following 
problem. 

PROBIF.M 1.1. Let D be a domain m the (x.y)- plane. and sup- 
pose F is a real-valued function defined on D . Suppose we 
are alven the values F^ = F(x^,y^) c>f F at some set of 
points (x^yp located in D, i = 1,2, ...,N. Find a function 
f defined on D which reasonably approximates F. 

This problem is, of course, precisely the problem of fit- 
ting a surface to given data. In many cases the domain D is 
a rectangle and the data points lie on a rectangular grid. 

There are, however, many practical problems (see the following 
section for some specific examples), where D is of unusual 
shape and where the data points are irregularly scattered 
throughout D. Thus, while we shall pay some attention to spe- 
cial methods for regularly spaced data, we are actually more 
interested in the general case. 

There are basically two approaches to handling Problem 1.1. 
First, we may try to construct a function f which interpolates 
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the data exactly; i.e., such that 
(l.l) f(x 1 ,y 1 ) F l , 1 = 1,2,. ..,N. 

This approach may be desirable when the function values at the 
data points arc known to high precision and where it is highly 
desirable that these values be preserved by the approximating 
function. 

The second approach involves constructing f which only 
approximately fits the data. This may be regarded as data 
smoothing and will be desirable when (as is often the case) 
the data are subject to inaccurate measurement or even errors. 
The question of whether interpolation or approximation should 
be used »’U1 not be discussed further here — this Is t problem 
which must be settled for the individual problem at hand. 

In discussing Problem l.l, it will be convenient to make 
a further distinction between those methods which are local in 
character (i.e., where the value of the constructed surface f 
at Che point (x,y) depends only on the data at relatively 
nearby points) and those methods which are global in nature. 
Thus, we discuss four categories of methods in sections 3-6: 

(1) global interpolation, (2) local interpolation, (3) global 
approximation, and (4) local approximation. In each of these 
sections we further subdi\ide the material according to the 
type of functions being used and the type of data (scattered or 
not) for which the method is suitable. 

In discussing methods which apply only to special arrange- 
ments of data points, we have two objectives in mind. First, 
the methods are of interest in their own right. More important- 
ly in terms of Problem 1.1, however, such methods can also be 
used in two- stage processes in which we first construct a sur- 
face g based on the scattered data, and then use g to gen- 
erate regular data for the construction of another (perhaps 
smoother or more convenient) surface f. Such two-stage methods 


will be discussed (along with several examples) in more detail 
in section 7. 

For many of the methods based on regular data and some of 
those for scattered data, error bounds are available to indi- 
cate how well smooth functions are approximated by the surface 
constructed. We do not have space to go into the extensive 
literature on error bounds. A simple test of how well a method 
will approximate smooth functions is, however, provided by its 
ability to reproduce polynomial surfaces exactly (that is', if 
F is a polynomial in x and y up to a certain degree, titan 
the surface f is identically equal to F) . For many of the 
methods we will be able to indicate the corresponding degree of 
exactness. 

In many of the applications of surface-fitting techniques 
(cf. the examples in section 2), the ultimate aim is to use the 
data to construct a contour map of the unknown function. Since 
F is known only at the data points, we must be content to con- 
struct a contour map for one of our fitted surfaces. In sec- 
tion 8 we discuss some approaches to accomplishing this numeri- 
cally. 

We close this introduction with a disclaimer--this survey 
does not include all possible methods for fitting surfaces to 
scattered data. For example, we have not discussed Fourier 
series methods, spatial filtering, and other such related sta- 
tistical techniques. In addition, the set of references for 
those methods which we have discussed are also not complete. 

My original intention was to compile as complete a bibliography 
as possible, but the sheer bulk of relevant papers and my in- 
ability to locate all of them convinced me to settle for less. 

I have opted to quote a fairly representative list of papers, 
including several other surveys. Further rcterences can be 
found by consulting these. 1 shall be very happy to receive 
information on references and methods I have overlooked. 


2, Examples 


In this section we shall quote several explicit examples 
of Problem 1.1 to emphasize the fact that unusually shaped re- 
gions and scattered data do arise frequently in practice. 

EXAMPLE 2.1. Petroleum exploration . In exploring for petro- 
leum, the contours of various underground layers of sandstone, 
shale, limestone, etc. can be important indicators of possible 
oil fields. Frequently, data on such layers is available from 
exploratory wells, which, however, have most likely been drilled 
at locations scattered randomly throughout some geographical re- 
gion of interest. To quote a specific example, Robinson, 
Charlesworth, and Ellis [166] consider precisely this problem 
for some data obtained from 7,500 wells drilled in Alberta. For 
another example of this type, see Whitten and Roelling [208], 
Problems similar to that mentioned in Example 2.1 arise 
frequently In cartography and submarine topography where the 
measurements represent actual elevations. In some cases the 
measurements must be taken from photographs or from sonar mea- 
surements and are usually subject to some measurement error 
(eg. see Kubik [125] for a discussion of photogramraetry) . 

EXAMPLE 2.2. Geological maps . There are a great many problems 
in Geology and the earth sciences in which the data arises from 
some other function of location besides actual elevations. For 
example, some geological variables of interest might include 
concentrations of various chemicals, specific gravity, electri- 
cal resistivity, grain size, texture, optical properties, iso- 
tope ratios, etc. To quote a specific example, Bhattacharyya 
[21, 22] discusses methods for fitting a surface to measurements 
(taken by airborne sensors) of magnetic potentials over a cer- 
tain portion of the Yukon. See also Bhattacharyya and Raychaud- 
huri [23] and Crain and Bhattacharyya [61], 
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The importance of 8urface- fitting methods in the earth sci- 
ences can be judged by the large number of papers in the area 
relating to various fitting methods. For a further list of 
problems and a discussion of some of the methods which have been 
applied, see the books of Bohrenberg and Ciese [31], Chorley 
[51], David [62], Harbsugh and Merriam [98], and Kerriam [140]. 
Recent survey papers include Whitten [203, 205] and Whitten and 
Koelling [207]. To add just a few more of the papers in the 
geological literature dealing with surface fitting to our list, 
we mention Anderson [7], Grant [91], Hessing, Lee, and Pierce 
[114], Holroyd and Bhattacharyya [115], Kubik [123, 125], Nor- 
cliffe [151], Reilly [162 ], Whitten [200, 201, 204], and Whit- 
ten and Koelling [206]. 

EXAMPLE 2.3. Heart potentials . In order to diagnose certain 
abnormal heart conditions, it is desired to make a series of 
several hundred contour maps of the heart potential field at 
time steps of 1/100 of a second throughout a heart beat. Data 
on these heart potentials can be obtained by fitting the patient 
with a shirt containing probes. Because of body geometry, when 
this shirt is flattened out it takes the nonrectangular form 
illustrated in Figure 1. Although the probes could be arranged 
fairly regularly in this domain, because of the added signifl- 



Figure l. Heart Potential Measurements 
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cance of frontal measurements, In practice more probes are 
fitted there than In the back. This example was brought to my 
attention by Ms. Patrlzla Ciarllni of Rome. 

Potential fields arise in many other applications. Wa 
have already mentioned Geology in Example 2.2. For some exam- 
ples in modelling plasmas see Buneman [40]. The problem arises 
in Biersack and Fink [24] in experimentally studying crystal 
structure using neutron bombardment. Data from waveform dis- 
tortion in electronic circuits can be found in Akims [5, 6]. 


3. Global interpolation methods 

In this section we outline several methods for solving the 
interpolation problem (l.l). 

3.1 Polynomial interpolation. (Scattered data) . The general 
theory of finite dimensional interpolation is, of course, very 
well known (e.g., see Davis [63]). Briefly, if [G^ er « N 
functions defined on the domain D, then the function 

N 

(3.1) f(x,y) Z a 0 (x,y) 

J=1 J J 

u 

will satisfy (1.1) if and only if [a.fj Is a solution of the 
linear system 
n 

(3.2) Za A (x ,y ) = F i = l,2,...,N. 

J--1 J J 1 1 1 


This system has a (unique) solution for arbitrary choices of 
data precisely when it is nonsingular. This depends on the 

choice of functions (0,}t and the location of the data points. 

J A N 

To illustrate this method, we may choose the {0^}” to be 

polynomials in x and y. Given N, there is some leeway in 

the choice of which powers of x and y 

with N = 3 


to use. For example, 
1, x, y or possibly 


one could use the functions 
2 2 

the functions 1, x , y , etc. Uhen S is of the form N - 


♦ I 


(d+l)(d+l), we might use the functions 


(0j(x,y) 3* = 



d ,d 
v=0 , u--0 


As simple as this sounds, there are some serious difficul- 
ties with polynomial Interpolation of scattered data. For open- 
ers, it is not so easy to guarantee that the system (3.2) is 
nonsingular. To give a very simple example, consider the case 
N < 3 with the functions 1, x, y. If the three data points 
happen to lie on a line, then (3.2) will in fact be singular. 
Even when (3.2) is nonsingular, it will often be the case (at 
least if N is moderately large) that the system will be ill- 
conditioned. Finally, as is well known, polynomials of even 
moderate degree exhibit a considerable oscillatory character, 
and the resulting surface (even though it is C°^ is often too 
undulating to be acceptable. The general problem of polynomial 
interpolation to scattered data is not usually treated in Nu- 
merical Analyses and Approximation Theory books (see, however, 
Kunz [126], Prenter [157], and Steffenson [186]). Some papers 
dealing with the question include Guenther [93], Thatcher [189], 
Thatcher and Milne [190], and Whaples [19 7}. Assuming the in- 
terpolant exists, error bounds have been studied in Clarlet and 
Raviart [52-55]. 


Let 


(3.3) 


_ r v p .m , n 

S> = span [x y 1 ’ _ 

m, n 3 ‘ V=0 , 


be the space of polynomials of degree m in x and of degree 
n in y. This linear space is of dimension (ra+1) (n+1) and 
is, in fact, the tensor product of the linear spaces P and 

t& 

^ . It is perhaps of interest to note that there always exists 
a (usually nonunique) polynomial p e } ^ which solves the 
interpolation problem (1.1), no matter how the data points are 
positioned, see Prenter [158]. 
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3.2 Polynomial Interpolation (gridded data) . We begin Chit 
subsection by defining what we mean by gridded data. Let 

(3.4) H = [a, bj a (c,d] 

be a rectangle, and let 


x 0 < X 1 < 

... < Vi = 

>0 < *1 < 

••• ' ’t* 1 ■ 


We suppose now that F is a function defined on H, and that 
we have the values of F at the corner points of the rectangu- 
lar grid defined by (3.5); i.e., 


(3.6) 


ij 


F(x i ,y j ). 


1 = 0,1,. . ,,k+l 
j = 0,1,... ,2+1. 


This Is a total of N - (k^2) (2+2) data points. 

It is quite easy to show that there exists a unique poly- 
nomial p in the class ^ ^ ^ (cf. the definition (3.3)) 
which interpolates the gridded data given in (3.4) -(3.6). In 
fact, p can be written down explicitly in terms of the one- 
dimensional Lagrange polynomials as 


k*l 2+1 

(3.7) p(x, y) Z Z F L (x)L (y), 
i=0 j=0 J 

where the (L^fx)}^ 1 and (y) )q + ^ at* the usual one- 

dimensional Lagrange polynomials associated with the interpola- 
k*l k+1 

tion points ( x i ^ 0 and ^ y i^0 ’ res P ecCtvel -y • Interpolation 

of gridded data by polynomials has been discussed in various 
books and papers--we do not bother with a long list. See e.g. 
Prenter [157] or Steffenson f 186 ] . More recently, there has 
been considerable work on Hermite and oeculatory interpolation 
in several variables; see e.g. Ahlin [3], Haussman [99,101,102], 
and Salzer [168-170]. 
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3.3 Shepard's method . In this subsection we discuss • method 
of Shepard (ISO) and some modifications of it. The method ap- 
plies to arbitrarily spaced data, and the interpolating function 
can be written down explicitly. 

Let o be some metric in the plane, for example the usual 
distance metric. Given a point (x,y), let r^ * p((x,y), (* l ,y i )l 
for i = 1,2, ...,N. Let 0 < m < “. Then Shepard's interpola- 
tion formula Is defined by 


(3.8) f (x, y) = 



z — 

r i M 


when 


r i 0, all 


when r. * 0. 


The formula (3.8) is defined for all points (x,y) In the 
2 

plane R . It is clear from the definition that it interpolates 


F t at the data points (x^y^, i = 1,2, ...,N. 


the values 

The value of f(x, y) at nondata points is obtained as a weight- 


ed average of all the data values, where the 1 


th 


measurement 
from the point 


is weighted according to the distance of (x,y) 

(*!.' y i> * 

Ue shall briefly recount some of the properties of Shep- 
ard's formula. First, by converting all of the terms to a 
common denominator, it can be shown that 
N 

( 3 . 9 ) f (x,y) = Z FA (x, y), 
i = l 1 1 

where 

N 

17 lr (x,y )] 4 

J-l J 

Ail 


(3.10) A t (x,y) 


N N ’ 

Z Z (r (x. >) 
k-1 £=1 e 
£tk 


i = 1,2, . . ,,N 


These functions satisfy 
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(3.11) A^(x^,y^) * i#j = 1>2, ...,H. 

The representation (3.9) is numerically more stable than the 
original formula (3.8). 

In view of its definition, ve see that the function f(x,y) 
constructed by Shepard is rot a simple polynomial or rational 
function. It is clear, however, that except for the points 
(x^,yi>, it is analytic everywhere in the plane. Its behavior 
in the vicinity of the data points (x^,y^) depends on the size 
of u. It can be shown that for 0 < u < 1, f has cusps at 
these points. For 1 < p, f has flat spots at the data points 
(i.e., the partial derivatives vanish there). We also observe 
the interesting property that 

(3.12) min F < f(x,y) < max F . 

lSi*N 1 ISiSN 1 

We may also note that if the data came from a constant function, 
i.e., F^ = c, i = 1,2, ...,N, then f is also the constant 
function f • c. 

We now comment on the choice of p. To get smooch surfaces 
without cusps, it is desirable to take 1 < p. On the other 
hand, if p is relatively large, then the surface tends to be- 
come very flat near the data points and consequently quite steep 
at points in between. Experiments (cf. Gordon and Wixom [90], 
Poeppelmeir [155], and Shepard [180]) seem to indicate that a 
choice of M = 2 is perhaps a good tradeoff. ([155] contains 
several examples showing the behavior as a function of p.) 

There are several drawbacks to Shepard's method (3.8), as 
pointed out by Shepard [180] Mmself. First, if N is large, 
then there is a very considerable amount of calculation in- 
volved in evaluating f(x, y) at a particular point. Secondly, 
the weights are assigned on the basis of the distance of points 
from (x,y) only, not their direction. Finally, the flat spots 



I 

I 

I 

I 

I 

i 

I 

I 

1 


f 

I 

I 

I 
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In the neighborhood of the data points is somewhat disturbing. 

The first of the-.e objections can be met by defining a local 

version of the formula, which we shall do in section 4.5. It is 

possible to construct an analogous formula which accounts for 

direction. For details, see Shepard [180]. Finally, we briefly 

discuss handling the flat spots. 

Suppose in addition to the function values F^ at each 

point (x^y^) we also have estimates FX^ and FY^ of 

F (x,,y,) and F (x,,y.). Then we may consider the function 
x i i y i i 

N 

(3.13) f (x, y) = ^A i (x,y)[F t + (x-x^FXj. + ty-y^FYj. 

It is easily checked that this fur_tion also interpolates, and 
that 

(3.14) f x (x t , yi ) = F* i f y (.x t , yi ) = FY 1 , i = 1,2,...,N. 

This property may be expressed in the assertion that if the 
data F^,FX^,FY^ came from a plane surface F, then f will 
exactly reproduce this surface. To use formula (3.13) in prac- 
tice on the data-flttlng Problem 1.1, we have to carry out a 
two-stege approximation process in which the first stage con- 
sists of some method for estimating the slope at each of the 
data points. 

It might be of practical interest in some cases to con- 
struct still a more sophisticated version of Shepard's formula 
which would exactly reproduce higher-order polynomial surfaces. 

One approach to doing this is to use the following lesaa. 

LFMMA 3.1. (Barnhill (15]). Let P and Q be linear projec- 
tions of some linear space of functions 9 into itself. Sup- 
pose that Q exactly reproduces the linear subspace E c 
t-g-. 


(3.15) Qp = p, all p e E. 
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In addition, suppose that { X^ is a set of linear func -* onals 
on 9 , and that 

(3.16) A Pf = X^f, all f c y, i = 

Then th~ Boolean sum projector • 

(3.17) P ©Q - P ♦ Q • PQ 

envoys the function precision of Q (i.e,, reproduces E) and 
the Interpolation properties of P (i.,e.. (3.16) also holds 
for P ©q) . 

This result permits tl.e construction of interpolation 
schemes using Shepard's formula which reproduce higher-order 
surfaces. For an example, see Poeppelmeir [ 15 5 J where Shepard's 
formula is combined with a certain local interpolation scheme 
which reproduces quadratic surfaces. In closing this section 
we note that Shepard's formula can also be inter reted as aris- 
ing from weighted least squarcs--see section 5.1. 

3.4 Spline interpolation (scattered data) . Suppose X is a 
linear space oi "smooth" functions defined on the domain D, 
and let 

(3.18) U (f ( X: f^x i ,y i ) - F , i = 1,2, 

U is the set of smootn functions which interpolate. Now sup- 
pose that 0 is a functional on X which measures the smooth- 
ness of an element in X--the smaLler 0(f) is, the smoother 
f is. Then we may consider the following minimization problem: 

(3.19) Find s c II such that 0(s) = inf 0(u). 

ut_U 

The function s will be the smoothest interpolant, and in view 
of the similarity with classical spline approximation, s is 
caileu a spline function interpolating F. The basic questions 
concerning spiine interpolation center around existence, unique- 
ness, characterization, and construction. A quite general 
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abstract theory of spline Interpolation has been built up (see 
eg. Laurent (12 7 J and references therein). In this section we 
quote some specific examples which can be used on Problem 1.1. 

Where X is a semi-Hilbert space, 0(f) = ||f||, where ||« j| 
is a seminoma on X, and N - (f c X: ||fj| = 0 }, it is possible 
to show (under some additional mild conditions on X, see Duchon 
[72, 73]) that problem (3.19) always has a solution which is 
unique up to an element in N. Moreover, it can be shown that 
there exists a reproducing kernel K defined on DxD such 
that 

*1 d 

(3.20) s(x,y) = L a K((x,y) ; (x ,y )) s E b p (x,y), 
i=l 1 LI 


where { ^ is a basis for N. Moreover, the coefficients 
(a^) and [b^] can be determined from the linear system of 
equations 


(3.21) 


N d 

E N((* j ,y i );(x i ,y 1 ))a 1 + ~ b i p i* x j ,y j ) =F j' J=1 » 




£ °» k " 1,2, ...,d. 

irl 1 * 1 1 


The development with semi-Hrlbert spaces in Duchon [72,73] 
is an extension of earlier work of Atteia [10-12] and Thomann 
[192-193] using Hilbert spaces. The essential difficulty in 
applying the general results is the construction of an appro- 
priate reproducing kernel. We turn now to some specific exam- 
ples. 

Suppose X is the space of all functions on the rectangle 

D = H (cf. (3.4)) which have (distributional) derivatives up to 

2 

order 2 which lie in L (10 . For f e X, let 
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The reproducing kernel in this case can be written down as an 

infinite series involving sin and cos, and the space N is 

spanned by 1, x, and y. Similarly, if we replace H by the 

unit disc UD, the kernel can be computed as an infinite series 

(see Atteia [10-12] and Thomann [192-193]). Thomann considers 

computation of these splines by approximating the infinite 

series--POKTRAN programs are also included. 

If we replace the bounded sets H or 'JD by the entire 
2 

plane R and introduce an appropriate space X, it is possible 

to obtain explicit expressions for t a r< Producing kernel. This 

is the content of tkichon [72, 73]. In paiticular, let H 

2 

be the set of all tempered distributions f on R whose 
Fourier transforms f satisfy / 1 f 1 1 S dt < ». Let X mS de- 
note the set of all iunctions which have derivatives up to 
order m lying in H S . Our first example concerns the space 
X 20 . If we choose 0 as in (3.22), then the interpolating 
spline solution of (3.19) is of thi form 
N , 

(3.23) s(x,y) = L a i r i (x,y) log (r t (x,y)) + b^ + b 2 y + bj, 

2 2 k 

where r JL (x,y) = [(x-x^) + (y-y t ) 1 . The coefficients are de- 
termined from the system (3.21) with d = 3, N = span [l,x,y), 
and K(z,w) = |z-w| log(z-w|. Duchon refers to this type of 
spline as a thin plate spline since the expression 6 relates 

to the energy in a thin plate forced to Interpolate the data. 

2 

This spline belongs to C(R ) . 

21 

As a second example, suppose we consider X - X . In 
this case the solution of (3.19) with 0 given by (3.22) has 
the form 

N 3 

(3. 24) s(x,y) ? a i (r i (x,y)) 4 - b 2 y 4 b^. 

Here K(z,w> =. |z-w | . Duchon [72, 73] refers to these splines 
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as pseudo-cubic splines because of the analogy with the cubic 
splines in one variable. They belong to C (R) . Pseudo quint lc 
splines etc. are also considered in Duchon [72,731. 

A similar program has been carried out by Mansfield [133- 
137) for some spaces of smooch functions defined on a rectangle 
H. In [136 ] she considers a space of functions T P>,n (Q;,P) , 
where m and n are positive integers and a < a < b, c<P<d. 
This space is actually defined by completion of a set of tensor 
product functions with respect to an appropriate inner-product, 
and we do not want to define it precisely here. A function 
f c T m ' n (a, P ) has the following properties, however: 
r 


(3.25) < 


<i,J) 

£ C (H) , i = 0 

A, 

. . . ,m- 1 

and j = 0, l, . . . ,n 

(s-J- 

l ’ ^ (x,:0 c AC [a 

,b) 

and f^ S 

J ’ j) (x,t3) cL 2 [a,bj ; 





j = 0,1, . . .,n-l 

(i,*- 

l_l) (Of, y) cAC[c 

,dl 

and f (1. 

s-i) (a,y) £L 2 [c,dl 





i = 0,1,.. . ,m-l 

(m-1, 

n_l) t AC(H) and 

c (“» n) 

t L 

l 2 (h), 

AC 

stands for the 

space ot absolutely continuous 


tions and where s - m-tn. By constructing an appropriate re- 
producing kernel, she is able to solve problem (3.19) with 


(3.26) 0(f) 


//If ( m,n ) ]2 i N f {t (5 " J '- i \x,^) J 2 dx 
H j-0 a 


m-l d , 

♦ / U U ’ ' («, y) Tdy. 

i 0 c 


In [1331, Mansfield carries out a similar analysis for a 
space of functions R n ’ n defined on the rectangle H. Here 
R ’ [a,b] x [c,d 1 , where L,[a,b] is the usual Sobolev 

space of tunctions with absolutely continuous derivatives up to 
order m-l, and with f ^ c L*(a,b). By constructing an 
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appropriate reproducing kernel, she now solves problem (3.19) 
with 

(3.27) 0(f) = //[f (n ' n) ] 2 + £ / (x,c) ] 2 dx 

H j=0 a 

m-1 d .. . - 

+ E / [f (i,n) (a,y)] 2 dy . 
i=0 c 

The solution turns out to be a piecewise polynomial of degree 
2m- 1 in x and of degree 2n-l in y. It is also in 
C“* * Z (H) . For the particular case of gridded data, it re- 

duces to the tensor product of one-variable splines (cf. the 
following section). Other more general definitions of 0 are 
also considered (with minor modifications on the one-dimension ■: 1 
integrals) . 

A more algebraic approach to constructing multidimensional 
spline functions (which also involves certain kernel functions) 
has been taken by Schaback [173-174}. His two-dimensional ker- 
nel function is obtained as a tensor product of one -dimensional 
kernels. 

3.5. Spline Interpolation (gridded data) . The problem of con- 
structing interpolating splines in two dimensions with gridded 
data as in (3. 4) -(3. 6) is, of course, a special case of the 
general problems discussed in subsection 3.4. The development 
of the gridded data case predated the more general development 
and, moreover, is considerably simpler. There are a great many 
papers on two-dimensional polynomial splines and generaliza- 
tions. We do not have space here to discuss all of them in de- 
tail. We shall be content to quote some of the papers and to 
give a somewhat more complete discussion of polynomial splines, 
which are the most widely used splines for this problem. 

Some early papers dealing with two-dimensional interpola- 
ting splines include Birkhoff and d~ Boor [26], Birkhoff and 
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Garabedian [27j, Price and Simonson [159], and Theilhcioer and 
Starkweather [191], In [26] certain bicubic splines were intro- 
duced which were later studied in detail in de Boor [32]. The 
problem was to minimize 

(3.28) / / [f (2 ' 2) (x,y)] 2 dxdy 
a c 

over all appropriately smooth functions on the rectangle H 

which interpolate the gridded data (3.4)- (3.6). It was found 

that the solution of this problem was a certain bicubic func- 

2 

tion with global smoothness C (H) . This problem was genera- 
lized to minimizing 
b d 

(3.29) 6(f) = / / [f* B » n, (x,y)rdxdy , m = 2p, n = 2q 

a c 

in A-hlberg, Nilson and Walsh [1,2], whose solution involves 

certain higher-order polynomial splines. Since they are widely 

used, we give a short algebraic treatment here. 

The points and (yj)g + * define a partition of the 

intervals [a,b] and [c,d], respectively (cf. (3.5>>. Suppose 

now that x, <...<x.<a<b<x, -<...< x_ , and 

1-m — — -1 — — k*2 — — x+ra-1 

y. <...<y,<c<d<y, »<...<y, .are chosen aibi- 
i-n “- 1 “* l + i+x \- 1 

trarlly. Let ^ t * le B-spiines associated with the 

x-partltion, and let the B-splines associated with the y-parti- 

tion be denoted by {N^(y))^ n - For a complete discussion of 

B-splines and their properties, see de Boor [36] in this volume 

(or [33]). Let 

(3.30) N (x,y) = N“(x)N^(y), i - l-m,...,k and J = 1-n 
The linear space 

(3.31) * * span j (x,y) ) t ^ J = f. n 

is clearly of dimension (k-m)(i’+n). We may now construct an 
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element in (3.31) which interpolates to the gridded data. 

Since there are only (k+2)(/+2) data points on the grid (cf. 
(3 .4) -(3. 6), it is clear that if we use <S to interpolate, we 
have 


(3.32) (k+m) (i+n) - (k+2) (i+2) = (k+2) (n- 2) +(J+2) (o-2) +(n-2) (m-2) 

free parameters. Thus, to uniquely define a spline, one must 
add additional conditions. Recall that m = 2p and n = 2q. Then 
we might add the extra conditions 

g( v ,°)( , y ) _ s (V,0) (st ,y ) =0, j =0,l,...,i+l 

(3.33) U 3 K+i j v = p, ... ,m-2 

s (° > n) ( x y Q ) = s (0 ' M) (x i ,y i+1 ) =0, i=0,i,...,k + l 

u = q, ...,u-2 


(3.34) 


9<V,4) <VV = S<V ' M) <V W - s(V ' M) (x k*l ,y 0 ) 

= s ( ,M> (* k+l »y /+l > =0 » V = p,...,m-i 

u — j • • • f n* i • 


These are called the natural boundary conditions, and it can be 
shown that the system of equations 

k i 

(3.35) Z I «<, N , SW# = F co» a = °,l, ...,k+l 


i= 1-m j-l-n lj U ° 


(3 = 0,1,..., /+1 


coupled with the conditions (3 .33) - (3.34) provides a nonsingular 
system of equations for the coefficients {a^}. This system 
has convenient bandedness properties if the equations are ar- 
ranged properly. The resulting SDline is precisely the solution 
of the minimisation problem (3.29). The boundary conditions 
(3 .33) - (3.34) are the natural ones associated with the problem 
(3.29). However, it is also possible to specify lower-order 
derivative information along the boundary and also obtain a 
nonsingular system of equations. The resulting spline, called 
Type I, can also be shown to satisfy an appropriate minimisation 


problem. However, for data- fitting purposes, to use the Inter- 
polate with boundary derivative data one would first have to 
perform a first-stage approximation to find estimates for the 
required derivatives. 

The best-known case of the above spline interpolation is 

the case o = n = 4, i.e. r bicubic spline interpolation. In 

this case the surface constructed is a piecewise bicubic with 
2 

global smoothness C (H) . The natural boundary conditions set 
second-derivative values to 0. Programs for computing natural 
bicubic interpolating splines can be found in the DiSL Library 
{11 7) in FORTRAN. FORTRAN programs for Type I bicubic splines 
can be found in Koelling and Whitten (121], where the required 
boundary information is assumed to be input. An ALGOL program 
for computing Type I bicubic splines in which boundary data are 
automatically computed by fitting one-dimensional splines ap- 
pears in Spath (183]. 

Bicubic spline interpolation has been widely applied. For 
some references in the Geology literature, see eg. Anderson [7], 
Bhattachnryya [22], Holroyd and Bhattacharyya (115], Koelling 
and Whitten (121], and Whitten and Koelling (2061. 

Problem (3.29) has been widely generalized in the spline 
literature. Instead of minimizing ordinary derivatives, one 
may introduce general linear operators, and instead of dealing 
with point evaluation functionals, more general linear function- 
als may be permitted. To list some (but by no means all) papers 
dealing with such generalizations, we mention Arthur (8,9], 
Birkhoff, Schultz and Varga [29], de Boor [34], Delvos (65,66], 
Delvos and Scheropp [68, 69], Delvos and Schlosser [70], Fisher 
and Jerome [78,79], Kaussmann [100], Haussmann and Munch [104], 
Munteanu [143,144], Nielson [148,150], Ritter [164,165], Sard 
[171,172], Schoenberg [176], Schultz [177,178], Spath [184,185], 
and Zavialov [209-212]. On L-shaped regions and other polygons 
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see Birkhoff (25 ] and Carlson and Hall [44-491. 

Ve close this section by mentioning another direction of 
generalization which has led to a considerable development, the 
idea of spline blending. These methods are useful for construc- 
tion of a surface which Interpolates not only function values 
at Isolated points but on the grid lines themselves; i.e., 

(3.36) f(x,yj) = F(x,yj) a < x < b and j = 0,l,...,/*l 

ffx^y) = F(x t ,y) c < y < d and 1 = 0, 1, . . .,k*l. 

To use such blending methods one must have F defined on the 
grid lines. Thus, the methods could be of value as second-stage 
processes. We do not have space to go into detail on spline- 
blended methods. We refer to the recent book of Barnhill and 
Riesenfeld [20] for a collection of papers on the subject and 
for further references. See also the papers of Gordon [84-87] 
and Gordon and Hall [88] . Recently, considerable effort has 
gone into showing that blending methods also arise as solutions 
of appropriate variational problems; see the papers of Delvos 
[65], Delvos and Rosters [66], and Delvos and Malinka [67]. 

4. Local interpolation methods 

The interpolation methods discussed in section 3 were glo- 
bal in nature — that is, the value f(x, y) of the constructed 
surface at any given point (x,y) in D depends on the values 
of all of the data points. i 111 s generally means that to euuipulc 
a representation for f one has to solve a fairly large system 
of equations, and to evaluate f(x,y) one generally has to 
carry out a considerable amount of arithmetic. In this section 
we shall consider local schemes where the surface depends only 
on nearby data points. Then the construction will usually lead 
to (a possibly large number) of sma ll systems of equations, and 
moreover, the evaluation of the surface at a given point will 
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usually involve very little computation. 

Many of the schemes mentioned in section 3 can be made lo- 
cal in nature by the following simple approach. Suppose that 

d 

the domain D is partitioned into subdomains: D = U D . We 

* 

then seek a surface in the form 


(4.1) f(x,y) = (f^Xjy), (x,y) e i = l,2,...,d. 

A* 

To construct each individual f^, we suppose that are do- 

mains containing D^, which contain only points which are "near" 
. Then we use the data (and only the data) in to con- 

struct f^. Usually, we can choose = D^. In most cases 
the most convenient choices for the subdomains are trian- 

gles and rectangles. We discuss these two cases first. 


4.1. Triangular subregions (scattered data) . Suppose that we 
are given data at points P = (x , y ), i = 1,2, ...,N scatter- 
ed throughout the plane, and let D be the convex hull of 
these points. It is more or less clear that by drawing lines 
from point to point we can construct a set of triangles with 
vertices at the which partition D. It is also cletr that 

given any set of points, this triangularization of D is not 
usually uniquely defined (see Figure 2 below for two different 
trlangularizations of the same region) . Moreover, as the fig- 
ure shows, some trlangularizations are superior to others in 
the sense that they exhibit fewer of the less desirable long 
thin triangles. 



Figure 2. Triangularization 
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The design of an algorithm to divide a region into accept- 
able triangles with vertices at prescribed points is not as 
easy as it sounds. Two algorithms in the literature which are 
designed to give good triangularizations can be found in Caven- 
dish [50] and in Lawson [128]. 

The simplest approach to defining a local interpolating 
surface is to construct f^(x, y) to be of the fora a^+ajXt 
a^y in each triangle. The data at the three corners of the 
triangle determine the coefficients for that piece of f (the 
corresponding system will be nonsingular provided the triangle 
is nondegenerate) . This procedure produces a piecewise linear 
surface which, in fact, will be globally continuous. This last 
property follows from the fact that along the sides of the tri- 
angle the functions reduce to straight lines joining the ver- 
tices. This method has been used by several authors for data 
fitting, e.g., Lawson [128] and Whitten [206]. For some con- 
touring routines based on this local interpolation scheme, see 
section 8. 

If we desire to Interpolate several sets of data defined 
on the same triangularization, it may be more convenient to 
compute Lagrangian functions rather than to compute the surface 


in each triangle separately. In particular, it is clear that 
we can construct functions [0j(x,y)}^ with the property 

(4.2) ^(x^yp = 6 , i,j = 1,2, ...,N. 


These functions car be constructed as pyramids in such a 


way that the function 
surrounding the point 


If j ““a w wiiaj wu v. uv wa Daugava 

(*j,yj) (see Figure 3). In terms of 


these Lagrangian functions, tuc interpolating surface is given 


by 

(4.3) f(x,y) 


N 

X F 0 (x, y) . 

M J J 
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Figure 3. A Lagrange Element 

The Lagrangian approach to local interpolation is very 
reminiscent of the finite element method in which the solution 
of an operator equation is sought in the form ot a linear com- 
bination of a set of functions (called elements) with the pro- 
perty (4.2). (See e.g., Prenter [157], Schultz [179], or 
Strang and Fix [188].) There is no need to restrict the ele- 
ments to be piecewise linear functions--we may use higher-order 
polynomials, rational functions, or even more complicated func- 
tions. In fact, if we are careful in the construction, we may 
be able to construct elements with small support but higher 
global smoothness. 

There are a great many papers in the finite-element litera- 
ture concerned with defining convenient smooth elements (La- 
grangian functions with small support) . To mention a few, see 
Samhill, Birkhoff, and Gordon [16], Barnhill and Gregory [17, 

18], Barnhill and Mansfield [19], Bitkhoff and Mansfield [28], 
Bramble and Zlamal [39], Goei [83], Hall [94], Mitchell [141], 
Mitchell and Phillips [142], Nicolatdis [146,147], Zenisek [213], 
Zienkowicz [214], and Zlamal [215-217]. The books on finite 
elements of Aziz [13], de Eoor [35], Strang and Fix [188], and 
Whiteman [198] should also be consulted. 
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The construction of elements with higher-order smoothness 
becomes increasingly difficult. For example, it is shown in 
hansfield [137] that to get an element with support on the tri- 
angles surrounding and with global continuity C*(D), it 

la necessary to use polynomials of degree 5 at least. (Hatters 
are somewhat simpler on regular trlangularizations, see subsec- 
tion 4.2 below.) 

We close this subsection by mentioning that it is also pos- 
sible to perform interpolation using elements based on triangles 
to data which also Involves derivatives, or in cnalogy with the 
blending methods, to data which includes function values along 
the edges of the triangles. (See e.g., Barnhill, Blrkhoff, and 
Cordon {16], or Barnhill and Gregory [17,18].) These methods 
are not directly applicable to the scattered data Problem 1.1, 
but may be useful as second-stage methods. 

4.2. Regular trlangularizations . When the data is distributed 
such that the region can be triangulated into a set of congru- 
ent triangles, then it is extremely ad/antageous to use the La- 
grange approach. In particular, in this case we can find an 
element 0 with value 1 at (0,0) such that all ether elements 
are translates of 0. In this case, f takes the form 
N 

(4.4) f(x,y) = £ F 0((x,y) - (x ,y )) . 

J=1 J J J 

We illustrate this with a couple of examples. Suppose that 
we are given data at points chosen from the collection 

(4.5) = ((i, J)) jL ^j r2 U ((i+fc, J+5)^ j eZ , Z = [integers } . 

These points lie on the comers of a triangular grid as shown 
in Figure 4. 

It is shown in Zwart [218, p. 673} thet there exists a func- 
tion 0 e C i (R^) which is 1 at the origin and 0 at all other 
points in 0, and nas support on the shaded region in Figure 4. 



n* r- r 


5J 




This function is constructed as a piecewise quadratic polyno- 
nl A similar element has been constructed by Powell p56] 
(th ;ure on page 267 of [156] should be rotated 45° to see 
t’ f . 


To give another example, suppose that we consider the set 
of points fl 2 which lie at the vertices of the grid defined by 
equilateral triangles shown in Figure 5. 



Figure 5. Another Regular Triangularization 
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It la shown in Fredrickson (SI) that there exists a function # 
which has value 1 at the origin and value 0 at all other points 
in flj. The function is in C^R 2 ), consists of piecewise 
quartlcs, and has support in the region shown in Figure 5. 
Fredrickson also constructs a piecewise cubic element with the 
same support but which is only C (R^) . For right triangles 
see Carlson and Hall {44]. 

4.3. Rectangular subregions . In this section we suppose that 
we have data given at points lying 'n a rectangular grid as in 
(3.4)-(3.6), and consider local interpolation methods. The 
sinnlest approach h're (cf. the triangularlsafion case) is to 
construct a separate bilinear function f(x,y) = a^ + a^x + 

*3 y + *4 Xy in each subrectan 8 le > H i j = ix i ,x i+l* x * y j’ y j+l*' 
using the four comer values to determine the coefficients. 

Since the bilinear patches reduce to linear functions on the 
grid lines, the global surface is C(R) . 

Several authors have considered constructing functions on 
each of the using higher-orcer polynomials. This requires 

additional information in addition to the four corner values. 

For example, if one seeks a bicubic 

3 3 

(4.6) f (x,y) = h 7. a x y , 
i-0 J-0 J 


there are 16 coefficients to determine. These could be deter- 
mined by the four comer values, plus the values of f , f , 

x y 

and f at each comer. To determine these, one must perform 
stxse first-stage process. Fcr some approaches to this, see 
Akims f5), Hessing, et al (114], and Shu, et al ( IS1 J . A FOR- 
TRAN program for Ak<ma's method can be found in (&). Nonpoly- 
norelal patches have also been considered; e.g., sec Birkhoff 
and Carebedisn (27). 

The Lagrange (firlte element) approach can al-o be used in 
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Che case of rectangular gridded data. T n particular, if ve can 
construct a f-mction satisfying (4.2) with local support, then 
the surface f given by (4.3) will interpolate and the method 
will be local in character. As before, the Lagrange approach 
is especially convenient if the grid is regular, l.e., if all 


H 


ij 


sub rectangles 
that the 
points lie in the set 


are congruent. To illustrate this, suppose 


are actually the unit squares; l.e., the data 


(4.7) Oj = {(1,3)) i,J e Z, Z = (integers). 

To get a quadratic C* element, we may simply rotate the ele- 
ment of Zwart [218] considered in the last section by 45 degrees 
(cf. Figure 4), or we may take the element of Powell [156]. 


4.4. Parametric representations . The methods discussed m the 
last section is concerned with data given on a rectangular grid. 
By using parametric representations, it is possible to construct 
similar local interpolating surfaces for data given at the cor- 
ners of any partition of D consisting of quadrilaterals. In 
this section we briefly describe how this might proceed. 

Suppose Q is a particular quadrilateral subregion of D 
of interest. In addition, suppose tnat x(s,t), y(s,t), and 
z(s,t) are functions defined on the unit square U = [0,1] x [0,1] 
with the properties that as (s,t) runs over the boundary of 
U, (x(s, t) ,y(8, t) ) runs over the boundary of the quadrilateral; 
the four corners of U correspond to the tour comers of Q; 
and z(s,t) takes on the desired data values at the four cor- 
ners of U. In this case, the triple (x(s, t) , y (s, t) ,z (s, t)) 
provides a parametric representation of a piece of surface de- 
fined over Q interpolating the data. 

The problem of constructing parametric representations of 
interpolating functions has been considered in a number of pa- 
pers. Several papers on these methods and a host of references 
can be found m tne book of Barnhill and Riesenfeld [20]; see 
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also the survey paper of Shu et al [181]. Such surfaces are 
sometimes called Coon's surfaces, cf. Coons [59], and arc of 
considerable Interest In the field of coraputer-alded geometric 
design. To mention Just a few of the actual papers, see Ahuja 
and Coons [4], Earnshaw and Youllle [74], Ferguson [77], Hayes 
[107], Hosaka [116], and Kangeron [132], 

There also has been seme effort directed towards construct- 
ing elements (Lagrange functions) associated with other less 
regular subsets of the plane We mention, for example, the 
work of Clarlet and Ravlart [55], Wachspress [194,195], and 
Zlaaal [217] In which elements are constructed for domains 
involving curved edges. 


4.5. Local Shepard methods . It Is possible to modify the meth- 
od discussed In subsection 3.3 to make it local. For example, 
following Shepard [189], suppose we fix 0 R and define 


r 


(4.8) \jr(r) - < 


1/r 


0 < r < i. 


H ({ - 1) 2 , R/3<r<R, 

0 , R < r 


This function is continuously differentiable and vanishes iden- 
tically for r < R. Now with r^ as in (3.8), we define 

% 


L F, [f(r >r 
1=1 1 1 
N 

(4.9) f (x, y) = t Z ftfr.)]* 
1 i=l 1 


when 


r t t 0, all i 


^ , when = 0. 

2 

Formula (4.9) is defined at alt (x,y) in the plane R . 
By definition it interpolates the values F^ at the data 
points (x^yp, i = 1, 2, ...,N. The values at non-data points 
are obtained as weighted averages of the data values F^, but 
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only those which lie at points within a distance of R of 
(x,y). Thus, the fornula is local. 

To use this method in practice it is necessary to choose 
a reasonable value for R. The aim is to find R so that for 
every (x,y) a reasonable number of data points will fall in 
the disk centered at (x, y) of radius R . It would also be 
possible to let R depend on (x,y), i.e., to use different 
value of R in different subregions of D. 

5. Global approximation 

As mentioned in the introduction, frequently the data does 
not warrant constructing an interpolating function (e.g., be- 
cause of errors) . In such cases it may be preferable to con- 
struct a surface which only approximates the data. In thio sec- 
tion we discuss sane global approximation mechods. 

5.1. Polynomial least squares . The general theory of discrete 
least-squares fitting is very well known. To briefly review, 
suppose that (0^)” arc n given functions on D. Define 

N | n |2 

(5.1) 0 (a) = Z Z a 0 (x ,y ) - F I , 

1-1 1 J=l J J 1 1 M 

where a - (a^,...,a ) x is any vector in R n . Then the prob- 
lem is to find a* such that 

(5.2) ®(a*) = min $(a) . 

a 

The corresponding function 
n 

(5.3) f(x,y) ^ Z a* 0.(x,y> 

j = l J J 

is called the discrete least-squares approximation of the cats 
( F ^ ^ . Usually one takes r. considerably smaller than N. In 
this section we briefly discuss least squares using polynomials. 
Before doing so. hewever, we make a few general remarks about 
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solving the general least-squares problem. 

There are several approaches to solving (5.2). Perhaps 


the neatest is the case where the 
respect to the inner-product 


{tfj}" are orthonorraal with 


N 

(5.4) (0,y) = >- 0(x ,y )* (x ,y ) . 

i-1 1 1 * 1 


Then the solution of (5.2) can be written down explicitly as 
n 


(5.5) f(x,y) = ^^0j(x,y). 


J 


A second very well-known approach to solving (5.2) is via 
the normal equations 


(5.6) 
where 

(5.7) 


A* A a - A*F , 

X 

F = (F, , ...,F.,) is the vector of data values, 
I N 


A “ 


(tf j (x 1 ,y i )) 


n , N 

J=M=i* 


and where 


In some cases the normal equations are a perfectly acceptable 
way to compute least-squares approximation, but in other cases 
the system (5.6) may be ill-conditioned (or even singular--cf . 
the following subsection for spline least squares). This ap- 
proach is also not convenient should side conditions be desired 
(e.g., by imposing actual interpolation at some of the values). 
For more cn the normal equations, see any book on Numerical 
Analysis. 

A more modem method of solving least-squares problems is 
to use general matrix methods. Specifically, consider the ob- 
servation equations 


(5.8) Aa = F. 



I 

I 


It can be shown that by applying a series of matrix transfor- 
mations to this system, one can o’, tain a set of equations giving 
the vector a*. For a complete description of methods of this 



\ 
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i 

type see Lawson and Hanson (129] or Stewart [187], Matrix t 

methods are quite amenable to the adding of side conditions and ' 

can also be designed to take account of rank-deficiency of the 

) 

matrix A (which corresponds to the case of singular normal 
equations) . 

Polynomial discrete least-squares fitting has been widely 
used for fitting surfaces to data, both scattered and regular. 

Several authors have developed algorithms for pr,/nomial dis- 
crete least-squares fitting of scattered data by constructing 
orthonormal polynomials (e.g. by Gran- Schmidt orthonormalisa- 
tion) . See, for example, Cadwell and Williams (42], Crain and 
Bhattacharyya [61], and Whitten (201,202]. The latter contains 
a FORTRAN program. 

When the data a.e more regularly distributed, polynomial 
least-squares fitting can often be simplified. For example, if 
the data lie on a grid as in (3.4) -(3. 6), then the desired or- 
thogonal polynomials are simply products of the one-dimensional 
orthogonal polynomials associated with the one-dimensional inner 
products corresponding to {xj}q + and (y^]g f * respectively; e.g., 
see Cadwell (41] or Clenshaw and Hayes [56], as well as the sur- 
vey papers of Hayes [105, 108, 109} . 

There are also special methods for handling data which arc 
not on a grid but instead lie on parallel straight lines. For 
example, Clenshaw and Hayes [56] have developed methods using 
expansions in terms of Tchcbycheff polynomials (although the 
method actually only produces an approximation to the least- 
squares fit rather than the actual minimum) . 

Polynomial least squares can also be interpreted as multi- 
dimensional regression us practiced by statisticians, e. ,, 
see Cffroyroson [75]. For example, if we arc try<ng to fit a 
function in the form 
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dx dy . . 
f(x,y) = Z Z a x y , 
i-0 j-0 1J 


then by defining new variables by 

V n V = 0,1, ...,dx 

*V(dy+l)+p ~ x ? f M = 0, 1, . . . , dy 

we can write 

d 

f(x,y) = Z b z d = dxdy + dx + dy, 

i-0 1 1 


and the problem becomes one of fitting a linear function in 
several variables. 

We close this section by observing that in some cases it 
may be desirable to consider weighted least squares. In parti- 
cular, if we have positive weights > 0, i = 1,2, ...,N, then 
we may replace C> in (5.1) by 
N . n 


I (a) - Z 

w 


' F i 


It is interesting to note that the interpolation formula 
of Shepard discussed in section 3.3 can be interpreted in terms 
of weighted least-squares fitting. In particular, fix (x,y) 
in D, and let r^(x,y) be the distance from (x,y) to the 
point (x^jV^) as before. Now set - r and consider 
least-squares approximation by a constant c, using these 
weights. Then one easily computes that the least-squares choice 
of c is 


c = 


N 

N 

Z w . 


Vi 


z 


This approach was adopted by Pelto, Elkins and Boyd [152] (as 
pointed out to mo by Chuck Duris). 
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5 . 2 . Discrete least-squares fitting by splines . As outlined 
in the previous subsection, discrete least squares can be car - 
ried out with any finite set of functions. It is not surpris- 
ing that a number of authors have tried using tensor product 
splines. See, e.g., Halliday, Wall, and Joyner [96], Hayes and 
Haliiday [110), Jordan [119], Hanson, Radbill, and Lawson [97], 
and Whiten [199]. Hayes and Halliday have developed both ALGOL 
and FORTRAN programs. It is, on the other hand, perhaps some- 
what surprising that least-squares fitting with splines can be 
somewhat problematical. We briefly discuss the method. 

Suppose that H - [a,b] x [c,d] is a rectangle containing 
the domain D of interest. Let (x^)q + ^ anc * P art ^” 

tions of [a, b] and [c,d], respectively, and let [N^]^*^ 

be the tensor product B-splines discussed in section 3.5. We 
consider discrete least-squares fitting using these (k-m) (f+n) 
B-splines. 

To explain hew difticulties can arise with spline least- 
square fitting, we observe that It is quite easy for the matrix 
A in the observational equations (5.8) to be tank-dcf icient. 

On a trivial level this can happen if for some B-spline N^, 
none of the data points lies in its support. This deficiency 
can, of course, be easily removed by dropping this particular 
B-spline from the set being used to approximate . But rank de- 
ficiency can also occur in more subtle ways because of the 
local support properties of the functions. This problem can be 
overcome with properly designed algorithms. See Hayes and Halli- 
day [110] for a careful discussion of spline least- squares fit- 
ting. Lawson and Hanson [129] Include a general discussion of 
hov to handle rank deficient lea3 t- squares problems. 

If we operate in terms of the normal equations, then It 
may v- 1 1 occur that the normal equations are in fact singular. 
This is again due to the local property of the B-splines com- 
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bined vith the discrete Inner-product. Even when it is not 
singular, the set of normal equations can be Ill-conditioned 
(even though it is a relatively sparse matrix with a kind of 
repeated band-structure). 

Discrete least squares can also be carried cut with vari- 
ous finite dimensional linear spaces of blended functions. For 
an extensive study of such methods, see the dissertation of 
Doty [71). 


5.3. Discrete and t approximation , 
discrete least squares, we oay consider 
approximation problem: Given functions 
we seek a* so that 


Instead of performing 
the following discrete 
(0j)” defined on D, 


(5.9) 0(a) 


n 

V 


i=l M 




*.i 


is minimized. Alternatively, we may minimize 


n 

(5.10) 0(a) = max j £ a 0 (x ,y ) - F |. 

ISiSN J-l J J 1 1 1 

These are the usual L and l best approximation problems. 

Both of these problems can easily be reformulated as linear 
programming problems for the determinations of the optical a* 
(cf. Rabinowicz [160,161] or Rosen [167]). Reasonable choices 
for the would be low-degree polynomials if D is small, 

or possibly spline functions. 

Discrete approximation methods of this type have had rela- 
tively little exposure in the literature. For some results 
using tensor product spline- in the i problem, see Rosen. 

C*J 

The optical a* was obtained there by using the standard sim- 
plex method on the associated dual linear programing problem. 

The £ problem can also be solved by using Rcmcz-type 
algorithms. For an algori the which performs generalized 
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rational approximation (and thus can a 'so be used for polynomi- 
al approximations) see Kaufman and Taylor [120]. Theoretical 
considerations for Tchebycheff approximation In several vari- 
ables can be found in Collctz [53] or Weinstein [196], for ex- 
ample. 

5. A. Spline smoothing (scattered dsta) . In this section we 
consider some minimization problems similar to those discussed 
in section 3.4, but where the class of admissible functions is 
not required to interpolate and where the functional to be mini- 
mized includes a term measuring how close the function comes to 
fitting the data. To be more specific, suppose X Is a linear 
space of "smooth" functions and that © Is a functional on X 
which measures the smoothness of an clement in X. Suppose in 
addition that E is a functional defined on X which measures 
how well a function fits the data. Then the spline-smoothing 
problem 13 the following: 

(5.11) Find s c X such that p(s) = inf p(u), 

ueX 

where 

(5.12) p(f) = ©(f) ♦ E (f) . 

The abstract theory of spline smoothing has been well 
developcd| see, c.g., the book of Laurent [127] and references 
therein. To illustrate the ideas, we briefly discuss a couple 
of examples. We suppose as in section 3.4 that X is a oeml- 
Hilbert space and that 0 is a semnora on X with H = 

[f c X: 0(f) = 0). We also suppose that X is actually a 
function space defined on a domain D, and that the point eval- 
uators at ((x^yj)}^’ are bounded linear functionals on X. 

We define 

»• 

(5.13) E(f) = p Z [£(x ,y ) -Fl 2 , 

i-l 1 1 1 

where p is a fixed positive constant. Then it can be shown 
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(cf. Duchon [72,73]) that the solution of Problem (5.11) la a 
spline which can be written in the form (3.20), where now the 
coefficients arc determined from the linear system 
N d 

£ K((x^,yj) ; (x 1 ,y i ))a i + L +a^/p = F^, 

(5 . 14) J = 1, 2, . . . , N, 

d 

T* = L - l,2,...,d. 


As in section 3.4, the application of this method depends 
on constructing a reproducing kernel K. If 9 is chosen as 
In (3.22), Atteia [10-12] and Thomann [192,193] considered 
spline smoothing for spaces of smooth functions on the rectangle 
and on the disc (the latter even contains ALGOL programs) . 

2 

Duchon [72,73] considers similar problems defined on D = R . 

A similar spline-smoothing problem has also been consider- 
ed by Pivorarova [154], where 9 i3 taken to be 

(5.15) 9(f) = / / [ D 2 f ] 2 + [D 2 f] 2 . 

x y 


See also Kubik [123], 


5.5. Smoothing splines (gridded data) . In 3ectio- 3.5 we con- 
sidered several minimization problems whose solutions led to 
interpolating polynomial splines (and generalizations) . In con- 
junction with the development of interpolating splines for 
gridded data, there was a concurrent development of smoothing 
splines. For example, instead of minimizing the integral 9 
in (3.29) over appropriate smooch interpolating functions, we 
may minimize instead o(f) - 9(f) + pE(£), where E is given 
by 

k + l 2 + 1 - 

(5.16) E(f) -- 1 £ [£(x ,y ) - F J . 

1-0 j-=0 1 J 1J 

For results in this direction, see e.g. Nielson [149, 150 J. For 
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9 given by (3.29), the smoothing splines are again polynomial 
splines. Again, more general linear differential operators and 
more general linear functionals can be considered. 


5.6. Continuous least square s. The method of continuous least 
squares i3 not directly suited to fitting surfaces to discrete 
data, but it can be of use as a second-stage process, so we 
briefly review it. We suppose now that F is a function de- 
fined on D which we wish to approximate, and that are 

given functions on D. We define 


(5.17) <f,g) = // f(x,y)g(x,y)dxdy, ||f|| 2 = (f, f> 

D 

and 

(5.18) C(a) = || E a.(J. - F|i 2 . 

j = i j J 


The problem is to find a* to minimize 0(a). The solution is 
given by solving the normal equations 

(5.19) Aa = r, 
where 

A = (^i^j))? Jal and * - K(J 1 »F)>...,<0 n ,F>J T . 

For reasonably nice approximating functions it is often 
possible to compute the normal matrix exactly. In practice, 
the difficulty lies in evaluating the right-hand sides. Gener- 
ally a quadrature formula is required for this. One advantage 
of the method would be that if several data-fitting problems 
are to be solved using the same set of approximating functions, 
one can do the work of Inverting the normal matrix just once 
and re-use the result as many times as desired. 

Reasonable choices for the approximating functions include 
polynomials, or better yet, tensor product B-spiinea as in 
(3.30). Here the singularity problems do not crop up for the 
splines because we are integrating instead of sunraing over 
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finitely many points. The normal matrix in this case has a kind 
of repeated band structure. The entries can be computed exactly, 
e.g., by Gaussian quadrature (cf. de Boor, Lvche and Schumaker 
[38]). Uniform best approximation by tensor products of splines 
has also been considered, e.g., see Sommer [182], 

6. Local approximation methods 

As pointed out at the beginning of section U, there are 
many advantages which accrue if one uses local methods rather 
than global ones. In this section we discuss some local approxi- 
mation schemes. 

6.1. Patch methods . As in the. case of interpolation, the sim- 
plest approach to obtaining local approximation methods is to 
partition the domain ard to define a surface (patch) on each 
subdomain separately. In particular, supjose that D = UfD,}^, 
where are disjoint subsets of D. Then we may seek f in 

the form 

(6.1) f(x,y) = [f t (x,y), (x,y) c D £ , 1 = l,2,...,d. 

To construct the patch y), we night use the data available 

in the subregion D^. In certain cases, however, it may well oc- 
cur that no data at all are available in the set D, . In this 
case we may choose a somewhat larger set cf points "near" 

D^, and use the data in to construct For any given 

method, it should be possible to make the choice of adaptive 

so that the size cf is kipt as small as possible consistent 

with the amount of data des'rcd fjr the construction of f^. 

The patch method outlined above can be used with any of the 
approximation methods discussed in section 5. For example, one 
might choose to use pol -rials (of low order), and to do dis- 
crete least-squares approximation. Or, one might use or 

f approximation or sc.c other convenient space (e.g. opline 3 ) 







instead of polynomials, The main point ie to keep the oite c2 
each individual patch problem (and thus the site o £ the corre- 
sponding system of equations) small. We nay have to solve a 
lot of systems of equations, but each will be small and fairly 
we 11- conditioned. 

To illustrate how the adaptive feature night be impl«aer.iJt.r , 
suppose that the domain D ot interest 1ms been enclosed la * 
rectangle H and that a partition of H ic defined by A - 

u(H ij 3 i=o'j o » wich u ij = lx i- x i+i J x fy j ,y j^i 1 fcr ,JOC 
(6.2) a - Xq < Xj^ <. . .< x v+l = b, c = y 0 <y t <... < y /+1 “ 

Now suppose that we want to do discrete least-squares fitting 

using a patch of the form f^^Xiy) - a + bx ♦ cy on H^. 

In -his case it would be reasonable to require that at least 

3 pieces ot data should be used to construct f^j. If 

does not contain 3 picce3 of dtta, wc expand to H < ^ Py 

adding all borderirg rectangles. If this does r.ot contain 

3 pieces, we again add all bordering rectangles, etc. Wc tl*s-t 

compute the discrete least-squares polynomial using the data lx 

H. ., but then wc use the resulting function orlv in H .. r'mc 
ij i) 

process may be repeated to define each rcqi ired patch. Thi* 
kind of adaptive algorithm is %’ery easy to program. 

In using patch methods to get local interpolation Dethoct, 
we concentrated on methods using data at comers oi triangle* 
or rectangles, end by choosing appropriate forms for the pacc-.es 
it was possible ro get tnc individual patches .o maten cogetc.tr 
to give a continuous global surface (or with n re oonhiaticatrs: 
patches, even C^(D) or higher). Here, however, where the in- 
dividual patches are determined by approximation, it is not 
very likely that the patches will match up. and the global sur- 
face will generally not even be continuous. Fcr most applica- 
tions, this is a serious drawback. However, as ve slrali see It 
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section 7, patch approximation methods can still be very useful 
as first-stage methods. 

6.2. Direct local methods . In this section we discuss some 
local methods in which the approximating surface is constructed 
directly from the data without solving any systems of equations. 
It will be convenient to pose a more general problem than pre- 
viously considered. 

Let 9 be a linear space of functions defined on D, and 
N 

suppose that (7^)^ arc linear functionals defined on 9. Let 

(0^)^ be a prescribed set of functions defined on D. Then we 
are interested in approximation scheres of the following fora: 

N 

(6.3) QF(x,y) -- Z 7 F0 (x,y) . 

i=l 1 

We can think of thia as a surface- fitting problem where the 
deta are given by F = /V F, i - 1,2, ...,K. Given the data, 
we can write the approximation down immediately. 

We also observe that if the have support on small sub- 
sets of D, and if each also has support on the same set, 

then the formula (6.3) is loca l. For example, if we take 7.^ 
to be point evaluation at the point (x^y^) and Ci^ (x, y) to 
be a function with support in a neighoornood of (x^,y^), then 
the approximation formula simply becomes 
N 

(6. *) QF (x, y) -= Z T 0 (x,y). 

i-1 1 1 

This is very reminiscent of the Lagrange form of interpolation 
(cf. (A. 3)), but unless the 0 are taken to satisfy (A. 2), 

QF will not in fact be an interpolmt. For this reason, for- 
mulae of the form (6.A) (or more generally (6.3)) are sometimes 
referred to as quasi- in terpol.-.nts . Local quasi-interpolants 
of the form (6.3) con be constructed 3inply by defining the 
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N 

functions (0^ witn local supports. If each of these is 
continuous (or smooth), then QF will also be. 

Although a host of quasi- intcrpolants can be constructed 
as outlined above, considerable care must be exercised in order 
to get methods which give good accuracy (when the original 
function F is smooth) . As observed earlier, this is directly 
related to making the method exact for polynomials, i.e., such 
that QP - P for all P in some class of polynomials. 

To construct methods of the form (6.3) which apply to 
scattered data, it is necessary to construct appropriate {0^}^. 
While a host of methods can be constructed this way, it is not 
so easy to choose the 0^ to make the method exact for poly- 
nomials (which, as we remarked earlier, is directly related to 
how well the method will approximate smooth functions F) . To 
get methods which do have a reasonable degree of exactness (and 

a correspondingly good error bound for smooth functions), it is 

N 

easier to first choose the (0.)., and then try to find suit- 
able While this generally rules out using point evalu- 

ators at scattered data, it is possible to construct methods 
based on point evaluators at appropriate points, and such meth- 
ods can be useful as second- stage approximations. 

To illustrate these ideas, we consider construction of 
local spline approximation methods following the general treat- 
ment in Lyche and Schumaker [131], Suppose D is enclosed in 

a rectangle H, and that H is partitioned into subrectangles 

k L 

by a grid as in (6.2). Suppose that (N .}, , ’ are the 

i ) 1= i-m, l-n 

tensor product B-spltnes associated with this partition (cf. 
(3.30)). We are now interested in approximation schemes of the 
form 

k £ 

(6.5) QF(x,y) =• !f, £ \ ,FN (x.y). 

i-l-tn j= l-n 

In particular, we are going to consider the question of 
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constructing fonnulac* of this type which are exact for the class 
of polynomials S , with some fixed 1 < v < tn ond 1 < u < n. 
This problem has a very simple algebraic solution if we decide 
to construct each . in the form 


(6.6) X 


ij 


v 

z 

v=l 


Z a,, X* X* 

ijvn ijv iju 


where the (X* ) V , and (X^. ) u , arc linear functionals 
ijv v=i iju n-1 

which apply to functions of x and y alone, respectively. It 

v 

can be shown (cf. [131]) that giver, any {X^j v } and (X* } sa- 

tisfying mild independence assumptions, there exist coefficients 

{ct, , } such that the formula (6.5) will be exact for S 

ijVp ' vu 

In fact, these coefficients can easily be explicitly computed. 

To give one example, suppose 
r 


(x. ,+...*x, ,) 

i^l i 4m— i 


(6.7) { 


(m-l) 

(n-1) 


i = 1-n, ...,k 
j = I— m,..«,r. 


Then wc obtain 


k i 

(6.8) QF(x,y) = Z Z F({ , *) )N .,(x,y), 

i- 1-n j- 1-n 1 J 1J 


a formula which exactly reproduces the bllineir polynomials S'. 
This is the multidimensional (tensor product) version of the 
Variation Diminishing method of Marsden and Schoenberg; it was 
studied in some detail in Munteanu and Schunakcr [165]. This 
formula is closely related to the Bczier-type surfaces construc- 
ted in Riescnfeld [163] (see also Gordon and Kiescnfeld [89]). 

We should observe that the way formula (6.5) now stands, It 
may involve information on F uhich is taken fiom data outside 
of the domain D. This situation ran be rectified as follows: 



69 


ORIGINAL ?>CZp 
OF POOR QUALITY 


(6.9) 0 = {(i,j): support H o not empty) . 

Then it can be shown (131) that the ->e^hod 

(6.10) QF(x,y ) = ZI \ 4 FN (x,y) 

(i,j)cn 

remains exact as long as all functions are restricted to D. 

To get higher-order methods, depending only on point eval- 
uations, we proceed as follows. Choose 


( 6 . 11 ) 


X i T ijv < X i-*in’ V " 

y j < T ij M < y j-m' 4 = l,2, *‘*' U ' 


for i = 1-m, ...,k and ) - 1-n, Then if we take \j v 

to be point evaluation at x* and A y to be point cvalu- 

iJV ija 

y 

ation at t, . , the coefficients in (6.6) are easily computed. 

1 la 

Hints on where the t's should be placed within the support 
of the B-nplines are siven by the error analysis in [131]. 

We close this section with some historical remarks on the 
development of local approximation schemes in two dimensions. 
Early papers include Babuska [14], de Boor and Fix [37], and 
Fix and Strang [80l . For some methods involving triangular 
partitions, see Fredrickson [32], Quasi-interpolants were 
constructed in de Boor and Fix [37] using point evaluation 
data, but including derivatives. We have followed Lyche and 
Schunaker [131] where general linear functionals are consider- 
ed, and where in particular, methods can be constructed using 
only point evaluation of the function. (Local integrals etc. 
would also be possible.) The papers [37] and [131] both con- 
tain extensive error bound analyses. It is striking that these 
local spline approximation methods give optimal order error 
bounds for smooth functions. 


* \ - 


*'1 



Many of the methods we have discussed in this paper are 
only applicable when the data are regularly spaced (and in fact, 
many surface-fitting methods require specification of derivative 
data as well an function values) . Such methods cannot be ap- 
plied directly to the scattered data-fitting Problem 1.1. On 
the other hand, some of the most convenient local interpolating 
and local approximating methods which do work for scattered 
data produce surfaces which are not globally smooth (or even 
continuous) . Thus, it seems natural to consider the possibility 
of constructing two-stage processes in which the first stage 
uses the scattered data to construct an approximation g, while 
the second stage uses g to generate data for constructing a 
surface f (with desirable properties, such as smoothness). 

Since it is quite clear how various methods discussed in 
the earlier sections might be put together to yield two-stage 
processes, it will suffice to mention Just a couple of examples 
here. 


7.1. Interpolation/intorpolatlcn . Suppose that we want to con- 
struct a piecewise bicubic surface based on data given on a 
rectangular grid as in (3.4) -(3. 6). In each subrectangle 
the 16 coefficients of the bicubic f (cf. (4.6)) would be de- 
termined by the values of f, f . f , and f at each of the 

' x / xy 

fou 1 - comers. Now since our original data-fitting problem only 

specifies the values of the function at the grid points, local 

interpolation cannot be carried out directly. However, we can 

use the data to provide estimates for the values of f , f , 

x y' 

and f _ at the grid points fi.e., we construct g interpolating 
xy 

the data); then we can use local bicubic interpolation as a 
second stage. The tender will have no difficulty in imagining 
ways to produce estimate -> lor these quantities. For sene meth- 
ods which appear m the Literal .re. 


I 

i 

! 

I 


see the papers of Akima 
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(5,6], Koellmg and Whitten (121], and Spath (183]. 

7.2. Approximation /Interpolation . Instead of making the first- 
stage process interpolation as in section 7.1, it would also be 
possible to use an approximating process. For example, one 
might use least-squares polynomial approximation to construct 

a patch surface and then use some convenient interpolation pro- 
cess as a second stage. For on example of this type, see Heas- 
ing ct al (114] . 

7.3. Approximation /approximation . This combination is parti- 
cularly convenient if we are not concerned about getting an in- 
terpolating function. Both stages can be made local. To give 
an example, recently I have constructed an algorithm for fitting 
surfaces to scattered data in which the first stage consists 

of polynomial least-squares patch approximation (with adaptive 
choice of data--sec section 6), and where the second stage con- 
sists of direct local tensor product spline approximation. Both 
stages are local, and the final surface is a tensor product 
spline. Since the second stage is a direct method, it i3 very 
cheap to apply. Experiments with real-life data (e.g. from 
heart potentials, potential fields, and geological nnps--see 
section 2) have produced very promising results. Details, in- 
cluding an analysis of error bounds, will appear elsewhere. T 
have also tried alternate versions where the patches are con- 
structed as low-order polynomials which are best approximations 
in the £, or £ sense (via linear programming) again with adap- 
tive choice of data. The results were very similar. Finally, 

I have also experimented with computing patch approximations, 
followed by continuous least-squares tensor-product spline ap- 
proximation. Again, the experiments were promising. 


8. Font curing 

As indicated in the introduction, rrequently the goal m 
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fitting a surface f to data is to construct a contour map 
which approximates Che contour map of the unknown surface F 
which produced the data. In this section we discuss some methods 
for constructing contour maps of a surface f. 

8.1. Piecewise linear functions on triangles . When the func- 
tion f to be contoured is a piecewise linear function defined 
on triangles (and globally continuous), locating contours re- 
duces essentially to a matter of good bookkeeping. Indeed, if 
H is the height of the contour of interest, then it is easily 
seen that for a given triangle T with vertices, A, B, and C, 

(8.1) the contour does not pass through T if H < mln(f(A), 
f (B) , f (C) ) or if H > max(f(A),f(B),f(C)) 

and 

(8.2) the contour intersects exactly two sides of T otherwise. 

If case (8.2) holds, it Is easy to determine which two sides 
are intersected and, moreover, by using inverse linear interpo- 
lation between vertex values, the points on these sides where 
the contour crosses can be determined. Specifically, if, for 
example, 

f (A) < H < f (B) , 

then the contour crosses the line from A to B at the point 
on the line which is a distance of 


(f(B)-f(A)) 


! B - A | 


from A. Given the points on two sides of a triangle where the 
contour line crosses, we can oow draw the contour line since it 
is simply a straight line between the points. An algorithm to 
carry out this procedure requires enunciating the triangles and 
vertices and sonc kind of effective search procedure. There 


are several available in the literature. For ALGOL programs, 
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acc Heap [111,112]. (An earlier paper of Heap and Pink [113] 
contains a similar FORTRAN program but only for regular triangu- 
larizations . ) Lawson [128] discusses a similar algorithm. The 
algorithms mentioned include two possible approaches: (1) one 

may start with a triangle where it is known the contour inter- 
sects, and trace this contour as far as it goes, or (2) one may 
simply draw the contour lines in all triangles which have them. 

8.2. Piecewise bilinear functions on rectangles . Suppose now 
that the function f to be contoured is a piecewise (continuous) 
function on a rectangle partitioned into subrcctangles by a grid. 
Since f is linear m x or y on the edges, it follows that 
we can again determine whether a contour line of height H 
crosses an edge by inverse linear interpolation. There is in 
this case, however, a serious difficulty which does not arise 
in the case of triangles. It may happen that the height H 
lies on three or even four sides cf the rectangle. In this 
case, it is possible that two different contour lines pas3 
through the rectangle, and it is not clear how to interconnect 
the points (see Figure 6). 



Put another way, if we are following a contour and enter a rec- 
tangle as shown above in Figure 6 on the bottom line, then it 
is not clear whether wc should now turn right or turn left. One 
approach to designing an algorithm in this case is to simply 
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always go right, say, even though this may In the end be wrong. 
(If it is, we have to start over with a coarser mesh.) This 
technique was Incorporated In an algorithm by Heap (111,112)-- 
the paper contains a FORTRAN prograta. (An earlier ALGOL pro- 
gram can be found In Heap and Pink (113). 

A second approach to handling the ambiguity problem is com- 
pute an approximation to the value of f at the center of the 
rectangle (e.g., by taking the average of the four-corner val- 
ues.) and then to triangulate the rectangle. This amounts to a 
second-stage approximation process, and the surface contoured 
is no leuger f itself but an approximation g. This idea was 
programmed in ALGOL in Heap and Pink (113) and in F0P.TRA11 in 
Heap (111,112). 

Once the set of points for a particular contour have been 
found, there are a variety of ways of drawing a contour line 
through these points. One possibility is to simply draw 
straight lines between each of the points. The actual contour 
lines arc expressions of the form y = (a+bx) /(c-dx) in each 
rectangle. These are generally not straight lines. Hence, if 
smoother contours are desired, one uay use any one of a number 
of methods for drawing a smooth curve through an ordered set of 
points in the plane. For example, the curve could be computed 
in parametric form using one-dimensional splines. Another pos- 
sibility would be to use the Bezier methods with either Bern- 
stein polynomials or with B-splines (cf. Gordon "’nd Riesenfeld 
(89) and Riesenfeld [163 ]), although in this case the curves 
will not exactly go through the points. For other algorithms 
see Harlow and Powell [138] or McConalogue [ 1 3 9 J . 

8.3. Piecewise Quadratics on triangles . Suppose now that f 
is a piecewise quadratic defined on a triangular partition. In 
this case a contour line at height H passing through a trian- 
gle must be described bv a conic section. Such a section can 
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be represented in parametric torn as 

x(t) « (b Q >b lC b 2 t 2 )/(b 3 +b 4 t ♦ b 5 t 2 ) 

y(t) - (b 6 + b ? t + b g t 2 )/(b 3 + b 4 t + b 5 t 2 ). 


see Powell (156 j. Powell has promised an algoritna based on 
this observation. 

We turn now to some methods for handling general functions 
f on arbitrary domains D. 

8.4. A simple line-printer method . The following simple-minded 

method can produce reasonable- looking contours without excessive 

computation, and without recourse to a plotter. Suppose H Is 

a rectangle enclosing the domain D, and that ve partition H 

a3 11 --- UH^ with a rectangular grid as in (6.2). Let HL < HU 

be given real numbers. Finally, suppose that t , Is some 

i j 

point In H^ where f can be evaluated (perhaps one of the 
comers or the center). Define 

r 


(8.3) 


'U 


, if f(t ) < HL 
, If f (t ) > HU 


v , if HL + (v-l)h < < HL - vh, l<v<8. 


for all 1 = 0,1,... ,k and J = 0,1,... ,1 (where hr (HU-UL) /8) . 
The (k+2) by (£+2) matrix C contains only integers, and if It 
is printed out without cither horizontal or vertical spacing, 
we obtain a reasonable- looking contour cap of the function. A 
typical example is included in Figure 7. The method can be 
refined by using an alpha-numeric array C and core than 10 
symbols. It can also be refined by using a printer with appro- 
priate horizontal spacing so that each symbol occupies a square 
rather than a rectangle (e.g., cf. bur.er.an [40]). 


8.5. Threading on a tectangular grid . As in section 8.4. 
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be represented in paramettic form as 

x(t) = ( b 0 +b 1 c +b 2 t 2 )/(b 3 -b^t < b 5 t 2 ) 

y(t) - (b^ + b^t 4 bgt 2 ) / (b^ + b^c 4 b,. t^) , 

see Powell [ 156 J . Powell has promised an algorithm based on 
this observation. 

We turn now to some methods for handling general functions 
f on arbitrary domains D. 


8.4. A simple line-printer method . The following simple-minded 
method can produce reasonable- looking contours without excessive 
computation, and without recourse to a plotter. Suppose H is 
a rectangle enclosing the domain D, and that we partition H 
as H - UK^ with a rectangular grid as in (6.2). Let HL < HU 
be given real numbers. Finally, suppose that t is some 

point in where f can be evaluated (perhaps one of the 

coraers or the center). Define 


(8.3) 


'U 


0 

, if 

£( V 

< HL 


9 

, if 

f(t i)> 

> HU 

-- 


if HL + (v-l)h < f(t 

ij> < HL 

+ vh, i < v < 8, 



for ail i = 0, 1, . . ,,k and j - 0,1,,..,! (where h = (HU-UL) /0) . 
The (k+2) by (i^-2) matrix C contains only integers, and if it 
is printed out without either horizontal or vertical spacing, 
we obtain a reasonable- looking contour nap of the function. A 
typical example is included in Figure 7. The method can be 
refined by using an alpha-numeric array C and more than 10 
symbols. It can also be refined by using a printer with appro- 
priate horizontal spacing so that each symbol occupies a square 
rather than a rectangle (c.g., cf. Bunenan [40]). 


8.5. Threading on a rectangular grid . As In section 8.4, 
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1?72<7223333344444S5555K6667777788P8PB8e387777 
272222222???33334445555666677777e888388S888777 
72111111112222333445555668777778S8388888688838 
2711111111112723344455566677778898888999999999 
?1111ooooo1111??334455586677788R8S999399995S39 
2111oooooooa11??33444556677783B399999999999999 
2111ocooooooo11??33445566777C89999999999999399 
711ocooooooao11723344556677B889999999999399999 
?11ocoooooooo11?23344556677B8B93999S9999999999 
21 loooooooool 112234455566778889999999999999999 
lllooooooool 1122334455666778883999999939999999 
111oooooo111172334455566777E88889999999999S999 
1111111111112233445556667778888889999999999999 
1111111111727334445566877778889888889998888383 
1111111227233344555666777788888888888888888888 
2222222723333445556667778888369839999999999393 
2277773333344455566777783888339399999339999399 
222333333444455 c 667777338888999999999999999999 
7333333444455556667778883339399399999999993999 
3333334444455556677778883839999999999999999999 
3333334444555566677778838338999999399999999999 
2333334444555566577777836388339993999999999199 
2333334444555566667777788888838399999993399999 
7233334444555566666777778888388388899799999399 
7233334444455556666677777788338883888999999999 
2233333444455555566666777771788888888999999999 
2723333444445555556665667777778383888899999999 
7723333344444555555566666777777788883899999993 
7273333344444455555566666677777778888889999939 
227233333444445555556666666777777R388885899999 
7277333334444445555566666667777777888883888899 

Figure. 7. A Simple Contour Hap (Heart Potential) 

suppose that D is imbedded in a rectangle H vhich has been 
partitioned by a rectangular grid as in (6.2). Assuming that 
f is continuous, it is still possible to decide which of the 
grid line3 a parti >lar contour of height H crosses by examin- 
ing the end-points of each such line. Since f is not generally 
linear along such a line, ue cannot determine exactly where the 
crossing point is by linear inverse interpolation. However, if 
we are willing to evaluate f a few times along this line, we 
can estimate the crossing point quite accurately by bisection. 
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for example. Once a sequence of points on a contour has been 
determined, we may thread a curve through the points just as In 
section 8.2. 

This methoo docs have one serious drawback, however, — 
jus t a3 with the method discussed in section 8.2--, if we are 
tracing a contour it may happen that after entering a triangle 
there is an ambiguity as to which of two points to use to exit 
the rectangle. One could opt for an ad hoc rule or try the 
second-stage approximation described in section 3.2. For at. 
example of how this method works, see Falconer £ 76 J (based on 
Lodwick and Whittle [130 J), where it is applied to a surface 
constructed by local weighted quadratic polynomial least-squares 
approximation. Since bisection is involved, one should realize 
that in drawing contours with this routing the surface f is 
going to be evaluated a great many times. 

8.6. Threading on a triangular grid . An obvious cure for the 
ambiguity discussed in section 8.5 for threading on a rectangu- 
lar grid is to use a triangular partition in the first place. 
Then the bisection method coupled with a threading routine leads 
immediately to a contouring routine for general surfaces f. 
Strangely enough, I have not been ade to find anywhere where 
this method lias been suggested. 

I have made no effort to track down all available papers 
on contouring. A few which 1 did find and have not yet men- 
tioned are Cottafawa and le Molt (60], Dayhoff [64], and Pelto 
et al (152], There are many others. 

In some cases it nay be desirable to have a more graphic 
picture of a surface than a contour nap can provide. Recently 
there has been considerable cffoit devoted to computer methods 
for displaying surfaces on a scope or with a plotter. For some 
examples of output and a discussion of methods, see e.g. the 
book by Barnhill and F.iesenfeld [20 i on computer-aided design. 
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l£ an actual 3-1) picture is desired instead of just a perspec- 
tive, it is even possible to produce holographs. 
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Abstract 


This paper describes an algorithm for constructing a smooth computable 
function, f, defined over the surface of a sphere and interpolating a set of n 
data values, u^, associated with n locations, p-, on the surface of the 
sphere. The interpolation function, f, will be continuous and have continuous 
first partial derivatives. The locations, p^, are not required to lie on 
any type of regular grid. 
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1. Introduction 


The problem of constructively defining a smooth surface that interpolates data 
defined at scattered points in the plane has been treated in different ways by 
a number of authors. For surveys of this work up to 1977 see Refs. (2) and (7). 

We consider here the analagous problem for data defined at scattered points 
over the surface of a sphere. When data is defined over only a portion of the 
surface of a sphere it may be most reasonable to map that portion of the 
spherical surface to a planar region, using a C* mapping function, and treat 
the problem by an algorithm designed for the planar domain problem. However 
when the data is scattered over the whole surface, and it is desired to 
produce a interpolation function defined over the entire surface, it 
seems necessary, or at least very desirable, to deal with the problem directly 
in the spherical setting. In particular, there is no C* function that will 
map the entire surface of a sphere to a bounded planar region. 

2. The problem 


Let S denote the surface of the unit sphere in 3-space. Given points 

p , i =1 , ..., n, the problem is to construct a computable function f, 

' 1 

defined and having C continuity over S, and satisfying the interpolation 
conditions 


f(p i ) = u for i=l , ...» n 
2.1 Relevant properties of functions on S 


A function of f defined on S is differentiable at a point p^ in S if and 
only if there exists a 3-vector g^ satisfying 


(1) f(p Q +dp) - (f(p Q ) + gjdp) 

lim = 0 

|dpt » 0 sdp B 


Pq + dp e S 
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Let Tq denote the tangent plane to the sphere at the point pg. Since the 
perturbed points p Q +dp in Eq. (1) are required to lie in S, the normalized 
perturbation vectors dp/ jdpi approach the plane Tq as g dpi approaches 
zero. It follows that if a vector gg satisfies Eq. (1) then so also does 
any vector of the form g Q +h where h is orthogonal to the tangent plane Tq, 
i.e. where h is a multiple of the vector from the origin to Pg. 

To resolve this nonuniqueness of vectors gg satisfying Eq. (1) we will 
standardize or, the shortest such vector. This vector is distinguished among 
vectors gg satisfying Eq. (1) by the property of being orthogonal to the 
position vector from the origin to Pg, or equivalantly by the property that 
the point Pg+g lies in the tangent plane Tq. We will call this vector 
g Q the gradient vector of f at p Q . 


Note that the fact that f has a restricted domain, namely S, is an essential 
part of this definition, for example if f is the restriction to S of some 
function f defined in an open neighborhood of 3-space containing p Q it is 
entirely possible that f may be differentiable at p Q and have a unique 
gradient vector g that is different from the (minimal length) gradient vector 
gg of f. In such a case however gg will be the orthogonal projection of g 
onto the 2-D subspace parallel to the tangent plane Tq. 

Let U be a region of S containing Pg and not extending more than v/2 radians 
away from pg in any direction. Let k be the one to one mapping of points of 
U to their orthogonal projections in Tg. Let Ug be the region in Tq to 
which U is mapped by k. Define the function fg on Ug by 

ys) = f(r 1 (s)) 

Note that the point Pg is in both the domains of f and fg. If f is 
differentiable at Pg with gradient vector gg then also fg is 
differentiable at Pg with gradient vector gg. We will make use of this 
local equivalence of f and fg later in deriving an algorithm for estimating 
the gradient of f from discrete data. 
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We will say a function defined on S is in the class C* if there is a 
continuous 3-0 vector-valued function g, defined on S, such that for each 
point Pg c S g(Pg) is orthogonal to the vector from the origin to Pg and 
satisfies the condition ascribed to g Q in Eq. (1). 

3. Major steps of the solution method 


The approach to be described has the same major steps as the method for the 
analagous planar problem given in Ref. (6). These steps are 

1. Build a triangular grid on S having the given points p. as vertices. 

2. Estimate the gradient vector g^ at each point p^ . 

3. To evaluate the interpolation function f at an arbitrary point p in S: 

(a) Look up p in the grid to find the triangle containing p. 

(b) Compute f(p) by an interpolation method using the given function 
values u. and the estimated gradient vectors g 1 at the three 
vertices of the enclosing triangle. 

3.1 Data structures 


In the algorithms to be described the points p i will br represented by their 
cartesian coordinates. It will be convenient in the following to let the same 
symbol denote either a point or the 3-D vector from the origin to the point. 

In particular, points in S are represented by vectors of unit euclidean length. 

Each triangle will have an index number and will be represented by a set of 
six pointers identifying the three adjacent triangles and the three vertex 
points. Tfns is exactly the same data structure as was used in Ref. (6). 

If triangle t has vertices whose indices are A, B, and C in counterclockwise 
order, and whose adjacent triangle indices are a, b, and c with triangle a 
opposite vertex A, b opposite B, and c opposite C, the six pointers 
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representing triangle f would be stored in one of the following three 
permutations: 

a, b, c, B, C, A 

b, c, a, C, A, B 

c, a, b. A, B, C 

All access to these pointers is done via three very short subroutines. Thus 
the actual storage mode for these pointers is "hidden" from the rest of the 
program. By appropriate programming of these three subroutines the pointers 
can be packed to save storage. 

The array storage requirements of this algorithm are thus 

3n locations for the vectors p^ , i=i, ..., n. 

n locations for the data values u., i =1 , ..., n. 

12n locations for the triangle pointers. This is based on 6 pointers 
per triangle and at most 2n-4 triangles. This storage requirement 
can easily be reduced by packing. 

3n locations for the gradient vectors g^ i=l, ..., n. 

n locations for a permutation vector used only while building the 

grid. This storage could be overlaid by the gradient vector array 
or could be eliminated entirely by minor changes in the program 
design. 

3.2 Determinantal tests and grid look-up 


Let Pj, p,, and p^ be 3-vectors havmj unit euclidean length. Let 
Det(pj, p^, p^) denote the determinant of the 3x3 matrix whose column 
vectors are p^, p,,, p^ in that order. 

If a = Det(pj, p,,, p^) £ 0 then no two of the vectors form an angle of 
zero or r, and the three vectors do not all lie m a single plane through the 
origin. In this case a proper spherical triangle can be formed by connecting 
each of tne three pairs of points bv '■he snorter arc of the great circle in S 
determined by that pair of points, inus each arc will have length less than n. 
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This triangle divides S into two regions. The smaller region is to be 
regarded as the interior of the triangle. If & > 0 an observer traversing the 
edges of the triangle with the interior of the triangle to the left will visit 
the vertices in the order p^, p^, o^. If a < 0 the ordering would be 
reversed. We will always order the vertices of triangles so that a > 0. 

Let p^, p^, pj, be vertices of a proper triangle t in S with a > 0. 

Regarding q as a variable 3-vector in S, note that the quantity 

s ! = Del (q,P 2 ,P 3 ) 

is proportional to the distance of q from the plane determined by the vectors 
P 2 and p^ with the sign of s^ being positive if q is on the same side of 
the p 2 pj plane as p^ and negative if q is on the opposite side. Thus a 
point q c S is inside the triangle t if and only if the three quantities 

= Det (q, p 2 , p 3 ) 
s 2 = Det (p lt q, p 3 ) 

S 3 = Det (pj, p 2 , q) 


are all nonnegative. 

Our algorithm for finding a triangle containing a given point q consists in 
computing the quantities s^, s 2 , s 3 for some triangle t and then either 
accepting t as the containing triangle if all s. > 0 or else moving to the 
neighboring triangle across the edge opposite vertex p^ if s^ is the first 
of the test quantities found to be negative. 

If there is no neighboring triangle across this edge the search stops, 
returning this information. Otherwise the search continues by computing the 
test quantities in the neighboring triangle. 

Rounding errors in :omputing a 3x3 determinant causing inconsistent sign 
determination could conceivably lead to cycling in the look-up process or to 
the construction of topologically impossible edges in the grid construction. 
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Consider for example four points p^ t p 4 that lie in order along an 

arc of a great circle, the arc having length less than *. The true 

mathematical value of the determinant of the 3x3 matrix formed using any three 
of these vectors is zero. 

Using finite precision coordinates 3nd finite precision floating point 
arithmetic these determinants will generally net be computed as zero. A 
nonzero result does not in itself cause a serious problem but the possibility 
of inconsistency in the evaluation of related determinants can. 

To illustrate the hazard suppose that with p^, ..., p 4 as above the 
computed value of Det(pj, p^, p^) is positive and (p^, p 2 , p 3 ) is 

accepted as a triangle in the grid. Then suppose p 4 is tested for inclusion 

in this triangle. It is possible that all of the determinants Det (p 4 , 

P 2 , p 3 ), Det (p lt p 4 , p 3 ), and Det(p 1> p 2> p 4 ) might evaluate 
nonnegative. This would lead to the erroneous conclusion that p 4 is 
contained in the triangle (p^, p^, p^) and various topologically 
incorrect edges would be constructed to incorporate p 4 into the grid. 

Using a tolerance c such that all results between -e and c are treated as zero 
does not solve the problem. We have had good luck using double precision 
evaluation of the determinants and strict zero tests. We have also had 
success with single precision determinant evaluation if we randomized the 
order in which the points p. were considered for inclusion in the grid. 

One way to assure consistency while sacrificing some accuracy would be to 
truncate all coordinate values to a small enough number of bits to permit the 
determinant evaluation to done exactly. For example, on a machine carrying 
fourteen hexadecimal digits of significance in a double precision number, one 
might round all coordinates to the 2~^ bit. The smallest nonzero bit that 
could occur in the product of th^ee such numbers would be the 2“^ bit. The 
coordinates do not exceed one in magnitude so the same is true of their 
products. These products and the sum of up to six such products can be held 
exactly in a normalized floating point number carring fourteen hexadecimal 
digits. Thus determinants of 3x3 matrices could be computed exactly. 



104 


3.3 Constructing the triangular grid 

The convex hull of a finite set of points in the plane is the smallest convex 
polygon containing the entire point set. We need an analagous notion, which 
we will call the spherical convex hull, for points on the surface S of the 
unit sphere. 

Let P be a finite set of points in S. If there is no plane that strictly 
separates the origin from P, we will say the whole surface S is the spherical 
convex hull of P. 

Alternatively if there is a plane strictly separating the origin from P let C 
be the smallest convex cone with its vertex at the origin and containing the 
set P. The intersection of C with S will be called the spherical convex hull 
of P. This region will lie strictly within some hemisphere of S. 

A triangular grid with n vertices and covering all of S will have 2n-4 
triangles. A grid that convers a spherically convex proper subset of S and 
has n vertices and b boundary edges will have 2n-b-2 triangles. Note that 2n 
can always be used as an upper bound on the number of triangles. 

Our method of constructing a triangular grid using a given finite point set P 

in S as vertices will be a sequential process that alters a grid covering the 
spherical convex hull of some set of k points of P to obtain a grid covering 
the spherical convex hull of these k points plus one more. 

Algorithms of this type can be divided into (at least ) three subtypes 

(a) First find the boundary points of the (spherical) convex hull of P and 
construct a triangular gni for these points. Then in the remaining 
sequential part of the algorithm each new point is known to lie in some 
triangle of the current grid. 

(b) Preprocess the points of P into an ordering that assures that each new 

point will be strictly outside the (spherical) convex hull of the 

preceeding points. 
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(c) Do no prepreprocessing and be prepared for each new point to be either 
inside or outside the (spherical) convex hull of the preceeding points. 

With subtypes (a) and (b) one is adding extra code and execution time for a 
preprocessing stage in the hope of permitting the subsequent sequential phase 
to be simpler and execute faster. We have at different times developed 
algorithms for the planar problem representing each of these subtypes. The 
algorithm of Ref. (6) is of subtype (b). My present inclination is to prefer 
subtype (c) as I think it permits the total program to be simpler and probably 
is not significantly slower if in fact it is any slower. More specifically it 
does not require storage for a separate data structure to keep track of 
boundary points as was the case in Ref. (6). 

Our approach then is to form one initial triangle and then loop through the 
remaining n-3 points adding one at a time and modifying the triangular grid at 
each stage to cover the spherical convex hull of all the points so far 
considereo. Each new point may be either inside or outside the grid so far 
constructed. 

In the class of problems for which this method is primarily intended, i.e. 
problems in which the data is scattered quite generally over all of S, a stage 
will be reached at which the spherical convex hull is all of S. Thereafter 
all additional points will necessarily lie inside the grid so far constructed 
since the grid will cover all of S. The user can cause this full coverage of S 
to happen early by arranging that the first four points to be processed are 
located such that the tetrahedron with these four points as vertices contains 
the origin as a strictly interior point. The triangular grid based on these 
four points will cover all of S. 

Initially the algorithm seeks three points with which to construct the first 
triangle. The first vector p^ is accepted unconditionally. The remaining 
vectors are scanned for the first one whose inner product with p^ lies 
between cos 179° and cos 1°, i.e. between -0.S9905 and 0.99985. Pointers 
are swapped to relabel this vector as p^. 
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The remaining vectors are scanned io find one whose determinant along with 
Pj and exceeds 0.001 in magnitude. Such a vector is relabeled as p^. 

The vectors p^ and are then swapped if necessary to assure that 
Det(p^, p^, p^) is positive. This completes the construction of tne 
first triangle. 

We may now assume a grid based on k-1 points has been constructed and the next 
point, p^, is to be introduced. A look-up is done using the method 
described in Sec. 3.2. This look-up either finds a triangle t containing 
P k , or else finds a triangle t such that p k is outside this triangle 
relative to a side of the triangle beyond which there is no adjacent triangle. 

In the first case, the single triangle t having vertex points p^, pg, 
will be replaced by three triangles hi .ng vertex points (p k , pg, pj.), 

(P a. P k » P(0> and (p A , p 3 , p k ) respectively. The algorithm then 
does a grid improvement phase to be described subsequently. 

In the second possible outcome of the look-up process, the point p k is 
strictly outside the spherical convex hull of the proceeding k-1 points, and 
in particular it is cutsiae an edge of triangle t that constitutes a portion 
of the boundary of the spherical convex hull. In this case one new triangle 
will be formed by connecting p k to the two ends of the edge of t that gave a 
negative s,. value in the look-up testing (See Sec. 3.2). 

The algorithm next scans the current grid boundary points in both directions 
from the hpw triangle and connects p k to all other boundary points that 
result in tne creation of proper spherical triangles (See Sec. 3.2). The 
algorithm then does grid improvement. 

3.3.1. Gr.d improvement 


When two adjacent spherical triangles form a strictly convex spherical 
quadrilateral there arises the possibility of replacing these two triangles by 
the two that occur when the quadrilateral is partitioned by its other 
diagonal. Gne must establish a criterion for choosing between <.he two 
possible disecuons of a quadrilateral. 
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This issue was discussed for the planar case in Ref. (6) where it was shown 
that three differently stated criteria were mathematically equivalent. In the 
spherical setting a fourth criterion with considerable intuitive appeal can be 
formulated and it is easily seen to be equivalent to the "circle test" of Ref. 
(5). 

Let Pj, p ? , p 3 , and p 4 be the vertices, in counterclockwise order, of 
a spherical quadrilateral in S. Assume all four of the potential triangles 
( p l P ? P 3 ), (P 2 P 3 P 4 ). (P 3 P 4 Pj). and ( p 4 P A P 2 ) 

would be proper spherical triangles. One choice would be to connect points 
Pj and forming triangles (pj p,, P 3 ) and (p^ p 4 Pj) while 
the otter choice would be to connect points p 2 and p 4 forming triangles 
n 3 p 4 ) and (p 4 pj p 2 ). 

Consider the 3-0 polyhedron underlying the spherical triangular grid. If the 
four points under consideraton are not coplanar then one choice will give 
underlying planar triangular faces ’’hat could be faces of a convex polyhedron 
and the other choice '.-ill not. This therefore is our new criterion: a 

preference to make the underlying 3-D polyhedron convex. 

Another way to describe this criterion is to consider the unique line L from 
the origin that intersects both of the lines p^p^ and PpP 4 - If Pj > 
p^, p 3< and p 4 are not coplanar the two lines will intersect L at two 
distinct points. We construct the one of these two lines that intersects L 
fur th crest from the origin. 

We implement this test by computing d = Det ( Pp-Pg » p-j— Pj , P 4 -Pj) 
and constructing che line p^p 4 if d > 0 ano constructing if 

d < 0. Either line can be used if d = 0. 

After a new point, say p^, is connected into the grid, each edge that is 
opposite tn some triangle is a candidate for swappinq. Thus if there 
is a triangle Pj { P 2 P 4 and an adjacent triangle p,P 3 P ( . the edge 
PpP, will be replaced by the cage p k p 3 if Det (p 2 -p k , P 3 -Pj , 

P 4 -Pj } is negative. When an edge is swapped the edqes opposite in 
the two newly formed triangles become candidates for swapping. 
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3.4. Estimation of gradient vectors 


We assume a triangular grid has been constructed in S covering the spherical 
convex hull of the points p^, .... p n and having the points p Jt .... 

P n as vertices. We also assume the data values u^ u n (See Sec. 2) 
are available. It is required to estimate a 3-0 gradient vector g^ at each 
point p. . See Sec. 2.1 for the characterization of gradient vectors for 
this problem. 

Let p^ be a point at which a gradient vector g. is to be estimated. Our 
general idea is to do a least squares quadratic fit to data near the point 
p^ and then use the gradient vector of this fitted quadratic polynomial as 
the gradient vector at p. . We use a six term quadratic polynomial in two 
variables forcing interpolation to the value u^ at p^. Thus we need at 
least five neighboring points, and prefer more than five to obtain a local 
smoothing effect on the gradient vector. 

Let Q denote the set of points to be used for tne fit. We first place all the 
immediate neighbors of p i into 0. If the number of immediate neighbors is 

from 6 through 16 and if the maxtnx for the least squares problem passes a 
conditioning test then this set Q is used for the fit. If the number of 
points exceeds 16, excess points are discarded. If the number is less than 6, 
more nearby points beyond the immediate neighbors of p i are introduced. If 
the matrix condition test is not passed, more points, up to 16, are added. If 
the condition test still fails with 16 points, the least squares system is 
damped to bias the solution toward small values of the coefficients of the 
three second order polynomial terms. 

The fitting is se' up in a 1 ocal coordinate system determined by p^ A 3x3 
rotation matrix R is determined that transforms the position vector of p i to 
the vector (0, 0, 1). Thus the "nc.-fh pole" of the rotated coordinate system 
is at p^. 
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The same coordinate transformation is applied to all vectors in the fitting 
set Q. Generally these transformed vectors, having some proximity to Pp 
will all lie in the "northern hemisphere" of the rotated coordinate system 
i.e. their z coordinates will be positive. If any transformed vector (x,y,z) 
has z < 0 we arbitrarily replace it by (x/s, y/s, 0) where s = Sqrt (x 2 + 
y 2 ). This last step is just an expedient to do something definite in a poor 
situation. Data must be very sparse or poorly distributed to result in any 
points of Q being in the "southern hemisphere" of the rotated coordinate 
system. 

We ignore the z coordinates of these transformed vectors, using only their x 
and y coordinates in the fitting. This can be interpreted as projecting the 
points Pj of Q orthogonally onto the plane T that is tangent to the sphere 
at the "north pole", i.e. at p. . The polynomial model for the fit is 

c l x + c ? y + c 3 x2 + c 4 xy + c 5 ^ = u-u i 

The coefficients Cp ..., c^ of this polynomial are determined by a least 
squares computation. The ?-vector (Cp c^} is the gradient vector at p i 
of the fitted polynomial relative to the xy coordinate system in the tangent 
plane T. Using the observations at the end of Sec. 2.1 we take the 3-vector 
(Cp Cp 0) to be the gradient vector at p 7 of the (as yet unknown) 
interpolating function defined over the surface of the sphere. The inverse of 
the rotation matrix R is then applied to (Cp c 2 , 0) to obtain the 
representation of the gradient vector g. | in the original coordinate system. 

3.5. Interpolation in a single triangle 

In the planar case described in Ref. (6) we preferred the 9-parameter 
Clough-Tocher cubic macroelement (Ref. 3) as our interpo 1 ation method 
primarily for the following two reasons: 

(a) It is more econo.ni .al to evaluate than any other C 1 interpolation 
method of which are aware. Beginning with the rectangular 
coordinates of q and of the vertices, and the function values and 2-D 
gradient vectors at the vertices, our evaluation of this mterpolant 
uses 55 multiplies, 65 adds, and 4 divides. 
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(b) The interpolant at any point is simply a cubic polynomial in the 

cartesian coordinates (or in the barycentric coordinates), and thus it 
is easy to derive and implement an evaluation of the gradient of the 
interpolated surface if this should be desired. 

Unfortunately, the Clough-Tocher method depends strongly on properties of 
polynomials in cartesian coordinates over a planar region and does not seem to 
generalize for use over a spherical triangle. 

We will describe two methods for C* interpolation over planar triangles that 
do generalize to sphetical triangles. These both represent the interpolant in 
the form 

(2) f(q) = Wjf 1 + w 2 f 2 + W 3 f 3 » w^ + w-> + w 3 = 1 

where the w-'s are nonnegative weight functions depending only on q and the 
locations of the vertices, and the f‘s are interpolants depending in general 
on q and all of the data associated with the triangle t, and satisfying some, 
but generally not all, of the conditions for continuity across triangle 
edges. A very helpful analysis of convex combination formulas of this type is 
given in Ref. (4). 

As with the Clough-Tocher interpolant the requirement for C* continuity 
across edges is approached by establishing values of the function and its 
gradient along an edge that depend only on data at the two ends of the edge. 

Values along an edge are computed by Hermite cubic interpolation and the 
tangential derivative at any point on an edge is computed as the derivative of 
this hermite cubic interpolation polynomial. The normal derivative at any 
point on an edge is computed by linear interpolation using the derivatives 
normal to the same edge at the two ends of the edge. For q on an edge of a 
triangle to let F(q) denote the value and G( q ) denote the gradient vector 
defined by these interpolation methods along the edge. 



3.5.1. 


Planar Method 1 
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For any point q in the triangle t let f in Eq. (2) be defined by Hermite 
cubic interpolation along the line througn q parallel to the edge opposite 
vertex p.. . This function f has been called the BBG mterpolant or BBG 
projector due to its use in Ref. (1). See also Ref. (2), pp. 92-101. 

Function and derivative values for this interpolation are derived from the 

edge functions F and G defined above. The function f.(q) defined in this 
. 1 1 
way is C over triangle t, and f and its gradient match F and G 

respectively on all edges except that the normal derivative of f. on the 

relative interior of the edge opposite p i will generally not be consistent 

with G. 

Corollary 2.5 of Ref. (4), adapted to our present notation, states that if all 

f/s in Eq. (2) match F on the entire boundary of t, and w^q) = 0 for any 

i and edge point q for which the qradient of f evaluated at q does not 
match G(q), then f of Eq. (2) matches F and the gradient of f matches G on the 
entire boundary of t. 

Thus letting f i be the BBG mterpolant, it will suffice to require that w l 
have the value zero on the edge opposite p and be nonzero elsewhere on the 
boundary of t. This is conveniently assured by lettinq w be the 
barycentric coordinate of q that has the value zero on the edge opposite p i 
and one at p . Thus Eq. (2) specializes to 

( 3 ) f(q) = b 1 f 1 (q) + b 2 f 2 (q) + *>3*3(9) 

where the b i are the barycentric coordinates of q relative to the triangle t 

and the f 's are BBG interpol ants , each requiring two linear interpolations 

and three Hermite cubic interpolations for its evaluation. 

3.5.2. Planar Methud 2 

For any point q in the triangle t let f in Eq. (2) be defined by Hermite 
cubic interpolation along a line from vertex p through q to the opposite 
edge. This mterpolant has been called a side-vertex or radial mterpolant. 
(See Ref. 2, 0 . 101). 
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The function matches F cn the entire boundary of t and its gradient 
matches G on the edge opposite p^ but its normal derivative is not 
consistent with on the relative interior of the two edges adjacent to 
p.j. Again using Corollary 2.5 of Ref. (4), it suffices to define w. of 
Eq. (2) to be zero on the relative interior of the two edges adjacent to p i 
and positive on the relative interior of the opposite edge. This is 
accomplished by setting 


w. 

i 


b i+l b i+2 


/(b i+ i b i+2 + b i+2 b i + b. b. +1 ) 


where the b. 's are barycentric coordinates of q and the subscripts are to be 
evaluated modulo 3 to one of the values 1, 2, or 3. 


The function w i defined in this way has non-removable singularities at 
vertices p i+ ^ and p 1+2 since it is one on the relative interior of edge 
^i+1^+2 anc * zer0 on * he re ^ at1ve ir >terior of the other two edges (and at 
vertex p^. For mathematical definiteness we may define w. ) to have the 
value zero at p i+ ^ and one at p 1+2 . The sum w, + w 2 + w 3 is then 
one throughout triangle t, as is required for Corollary 2.5. In a computer 
implementation one would treat interpolation at a vertex as a special 
(trivial) case anyway, so the particular choice of definition of w.. at the 
P 1 + j and has no bearing on implementations. 

Thus Eq. (2) specializes to 

(4) 

! ( b 3 b 2 f l + b 1 b 3 f 2 + b 2 b 1 f 3 )/(b 3 b 2 + b^ + b^), for q ^ p^ p 2 , or p 3 
u ,, Tor q = p i 


where the b^'s are the barycentric noordinates of q and the f 's are 
side-vertex interpolators. Each f requires one linear interpolation and 
two Hermite cubic interpolations for its evaluation. 


3.5.3. 
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Generalization of Planar Methods 1 and 2 for spheri-al triangles 


The key in generalizing these two planar methods for use with a grid of 
spherical triangles on the surface S of the unit sphere is to replace all of 
the linear and Hermite cubic interpolations along line segments by the same 
type of interpolations along arcs of great circles in S. 

Let t denote a proper spherical triangle with vertex position vectors p^, 

P 2 , and p 3 , and let q be a point of S contained in t. Let t' denote the 
underlying planar triangle having the same vertices as t, and let q* be the 
central projection of q into the plane of triangle t', i.e. q' is the point i 
the plane of t' intersected by the line from the center of the sphere to q. 

When the look-up procedure of Sec. 3.2 finds that a given point q in S is in 
triangle t, it also retu-ns the three nonnegative numbers s^, s 2 , and 

We call these numbers unnormalized barycentric coordinates since the 
(normalized) barycentric coordinates of q' relative to the planar triangle t' 
can be computed as 


b, = s 1 /( s x + s 2 + s 3 ), i = 1, 2, 3. 

The intersection points between certain lines through q' and edges of t' 
needed for either of the two planar interpolation methods are easily 
represented in terms of the b 1 's and p 1 's. Thus the intersection betv/een 
edge with the line through q' parallel to edge P 1 +iP 1+ 2 bas 

position vector b 1 p i + (l-b^p^j while the intersection between edge 
p i + l p i + 2 with the line from vertex o 1 through q has the position vector 

^ b i+l p i+l + b i+2 p i+2^ b i+l + b i+2^ * 

These intersection points can then be centrally projected to S by normalizing 
their position vectors to have unit euclidean length. All of the linear and 
cubic interpolations called for in the planar methods are then done wth 
respect to arc length along great circle arcs in S obtained by central 
projection of the corresponding line segments in the planar triangle t'. 
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Recall that gradient data at each vertex p^, p^, and p 3 is represented 
as a 3-vector that is orthogonal to the position vector of the vertex. 
Gradient information generated at auxiliary points in either interpolation 
method is also represented a 3-vector orthogonal to the associated position 
vector. 

The verification that each of these two spherical triangle interpolation 
methods defines a C* function over 3 can be carried out in the same way the 
C 1 property of the planar methods is proved, for example in Ref. (4). Thus 
one observes that the function value and gradient vector at any edge point of 
a spherical triangle is determined only by data at the ends of the edge and 
thus will be consistent in neighboring triangles. The partial interpolation 
functions f. have the correct values at all edge points and gradient values 
that are correct on certain edges and wrong on others. The convex combination 
formula (3) cr (4) properly zeros out the functions where their gradient 
values are wrong and thus gives a function having the required boundary values 
and boundary gradients. 

4. Software implementing these algorithms 

Subroutines were written for these algorithms in 1979 using the JPL SFTRAN3 
structured Fortran language which is preprccessed Federal (ANSI) Standard 
Fortran 77. 

The time for grid construction for n points was orcportional to n A *“ for 
test cases in the range from 25 to 500 points. The RMS error in test cases 
using simple mathematical functions to generate data over relatively uniform 
triangular grids of various densities was proportional to h -5 * 4 * in test cases 
having maximum edge length in the grid ranging from 63° down to 9°. 

A count of the number of arithmetic operations required to do a single 
interpolation in a triangle gives t^e figures listed in Table 1. The planar 
Clough-Tocher method is included for comparison. For all methods the 
computation starts with cartesian coordinates for q, p^, p^, and p 3 and 
function values and gradient vectors at p^ p^, and p 3 . The weights used 
to combine the counts are arbitrary but plausible. They are normalized to 
cause an add plus a multiply to sum to one for consistency with operation 
counts measured in ,, nops B . 
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For Method 1, tests of the continuity of the interpolated function and its 
first partial derivatives across edges of the grid were made in two ways. 
Interpolated values and their first and second differences were computed at a 
sequence of equispaced points along a smooth arc. Paths were chosen crossing 
edges at various angles and crossing a vertex. These tests indicated 
continuity of the values and first differences with discontinuity of the 
second difference at edges. 

The other test of C* continuity involved reprogramming all of the code for 
Method 1 using a "U-arithmetic" package developed at JPL in 1971 based on the 
ideas of Ref. (8). (This is like the method of Ref. (5) without the benefit 
of a preprocessor.) In this approach the program computes a 3-0 gradient 
vector and a 3x3 Hessian matrix for every intermediate quantity and thus also 
for the final interpolated value. All derivative computations use 
matnematical ly correct formulas, i.e. not differencing. 

We found it necessary to reorder some computations to avoid severe artificial 
numerical instabilities in the derivative computations. After this reordering 
the results were consistent with continuity. 

We did not try a U-arithmetic version of Method 2. I would expect severe 
difficulties with this since the singularities of the w- *s at certain 
vertices (See Sec. 3.5.2) imply that some first partial derivations of the 
w^s can be arbitrarily large in a small neighborhood of a vertex. 
Mathematically these cancel out but numerically there would be large rounding 
errors. 

5. An application 


In Fetruary, 1932, this software was used at JPL in the study of gravity 
variation over the surface of the planet Venus. Data was available at many, 
but not all points, of a rectangular longitude-latitude grid. The missing 
data occurred in irregularly shaped regions determined by geometrical 
constraints of the observation and communication instruments. 
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Table 1. Operation counts for a single 
interpolation in a triangle 



Clough-Tocher 
Planar Method 

Spherical 
Method 1 

Spherical 
Method 2 

Factors 

for 

weighted 

total 

Add/Subtract 

65 

371 

352 

0.4 

Multiply 

55 

699 

450 

0.6 

Divide 

4 

81 

57 

1.2 

Sqrt 


24 

15 

3.0 

At an 


18 

12 

5.0 

Weighted 

Total 

(Flops) 

63.8 

827.0 

584.2 
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Using 2450 points at which data was present our program built a spherical 
triangular grid consisting of 4896 triangles. Missing data in the rectangular 
grid was then filled in by interpolation in the triangular grid. 

In the course of this work the scientists gained new insights regarding their 
data and we found and repaired a weak spot in our program. See the discussion 
of determinant evaluation in Sec. 3.2. 

6. Conclusions and remarks 


The efficiency of the grid building procedure, execution time in test cases 
being observed to be proportional to n^*^, is quite satisfactory. 

C 1 interpolation in a spherical trianole requires 9 to 13 more Flops than 
C* interpolation in a planar triangle. Modifications giving small 
reductions in the operation counts are known but it would be interesting if an 
entirely different approach could be found that might be more intrinsically 
related to the topology of the spherical surface and require significantly 
fewer Flops. 

Method 1 is more time-consuming than Method 2 by a factor of about 3 to 2 
since Method 1 uses nine cubic interpolations along arcs compared with six for 
Metnod 2. Analytic computation of gradients for interpolated values would 
probably be more stable using Method 1 than Method 2 because of the 
singularities in the w^ 's of Method 2. It would be interesting to make 
visual comparisons of surfaces generated by these two methods but we have not 
had the resources to make such comparisons. 

The programs appear to be robust and reliable. The use of the SFTR AN3 
structured Fortran language has been extremely helpful in keeping the code 
understandable. 

It should be noted that the use of the surface of a sphere as the domain is 
just a mathematical construct for dealing with the set of all directions in 
3-space from an origin point. Thus the methods of this paper are applicable 
to the representation of any bounded two-dimensional C* surface in 3-space 
that is "starlike" in the sense that there is some origir point from which a 
ray in any airection intersects the surface in at most one point and the ray 
is not tangent to the surface at that point. 
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Other two dimensional manifolds besides the plane and the spherical surface 
that may deserve investigation for scattered data interpolation include the 
surface of a cylinder or a torus. On a cylinder one may wish to admit 
triangles having two vertices at the same data point while on the torus one 
may admit triangles having all three vertices at the same data point! 
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SURFACES : R epresentation and Approximation 

The general: Surface Application 

"ree-form Sculptured 


Design of Surfaces 
Representation of Surfaces 


Examples: 

Lockheed 

Energy 


Cri teri a/characten sties: 

Fits given infornation/applicaticn 

Smoothness 

Shape fidelity 

Parametric representation 

Local »■'. global schemes 

Interacti ve design 

Interactive viewing 
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Outline cf 3D Surfaces 
Interpolation 

I 

Patches v. Points 

I 

Coons v. Bezier 

1 

(1) □ patches (2) A patches 

I 

Transfinite 

Boolean sums/correction surfaces 
Lofting interpolants 
Compatibility conditions 

Discretization -*■ Point methods 
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Approximation 

I 

Least Squares 
(J. G. Hayes, NFL) 


Shepard's Formula 


Preprocessors : 

1. Triangulation 

1st pass 
Optimization 

2. Gradient specification 


Triangular Interpolants 

Barycentric coordinates 
BBG 


Radial Nielson 
Symmetric Gregory 
Convex combinations 


(1) □ Coons' Patches 
First Problem : Interpolate 

to the 4 curves. 


ORIG’N'W- w 

OP pOCii QUALITY 


S. A. Coons 1964 


F "primitive" 

F general coordinate 
X(u,v) - 

y(u,v) parametric surface 
z(u,v)_ 


F(l,v) 




Solution to the 1st problem: Lofting interpclant 
F X F 5 (l-u)F(O.v) + uF(l.v) 

Univariate linear interpolant: 

f = f(u) - P x f - (l-u)f 0 + uf 1 . 

Bivariate F = F(u,v) ~ P^F •= the above. 


(O.fQ) 



Error F - P^F 
Idea- Mated. F-PjF 
and add this to P^F 


Solution: P£[F - P^F] 

does the job. 
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ORlGU'IAi- FAcI & 
OF POOR QUALITY 


Geometry vs. Algebra 


g - g(v) - P 2 g - (l-v)g Q + vg x 

PjCF-PjF) - (1-v) (F - PjF) (u,0) + vCF-P^Cu.l) 

~ (l-v)F(u,0) - (1-v) [ ( 1-u) F (0 , 0) + uF(l,0)] 

+ vF(u,l) - v[ (l-u)F(0,l) + uF(l,l)} . 

Overall approximation PF = P^F + P 2 F - P^P 2 ^ 

= (l-u)F(0,v) + uF(l,v) + (l-v)F(u.O) + vF(u,l) 

- [(l-u)(l-v)F(0,0) + u(l-v)F(l,0) + (l-u)vF(O.l) + uvF(l.l)] • 

Check interpolation: (PF)(u,0) = F(u,0) etc. 

transfinite interpolant, blending functions 
W. J. Gordon 1969 


C“ . 


PF = (P^ P 2 )F Boole, l sun 


w 


tensor product 


Bilinear ly blended Coons' patch 
Piecewise method - C® 

Practical applications: C^ or 




bicubically blended Coons’ patch. 



Ottawa. « 

OF POOR QUALITY 
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f =■ f(u) - P,f = hg(u)f(0) + h,(u)f(l) + ^j(u)f'(0) + K,(u)f'(l) 


F a F(u,v) 



(h n (u)h, (u)K n (u)K, (u) ]B 


V v) ' 

h,(v) 

F o (v) 

K,(v) 


vfoere 


B - rF(O.O) F(0,1) 
F(1,0) F(l.l) 

Fi oCO.O) F uo ( 0 . 1 ) 

r i.o (l ' 0) 


F 0il (0,0) F 0tl <0.1) 

F o,i a * 0) F 0.1 (1 ' l> 

F l,l (0,0) Fl.lW.W 

Fi.ia.Q) f, ,(1,1) 


F i F o.: 


F i,o! f i,i 


Twist trouble ("F, ," = 


Positions i Tangents 

I - 

Tangents J Twists 


Gregory's Square 


u - 


zh- 


3v3u 


( 0 , 0 ) + v 


3 2 T 
3u3v 


( 0 , 0 ) 


u + v 


Discretization •* Point ISethods 

F(u, 0 ) = F(u, 0 ) h h 0 (u)r(u, 0 ) +h 1 (u)F(l, 0 ) + Q ( 0 , 0 ) + V u)F l,0 (1 > 0) 
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ORIGINAL PASH IS 

OF POOR QUALITY 


Triangular Patches 


Rectangular detrains vs. non -rectangular domains 


Preprocessors: a. Triangulatlcn 

Algorithm: (1) Enforce given boundary 

Default: convex hull. 

(2) A triangulation ... fast. 

(3) Optimize: min max angle 

T T6T 

where T is the set of tri angulations . 


b. Gradient Specification 
Surface Design: use tangent handles 

Surface Representation: use triangular Shepard's Method (Little) or 

inverse-distance-weighted least squares (Franke) . 



Linear interpolant = bjF(Vj) + b 2 F(V 2 ) + b^FCV^) 
(Finite elements) 


Problem: 


Find triangular interpolants 
C 1 triangular Coons' patches 

Barnhill, Birkhoff, Gordon, 1969-73 


ORlCIX " L 

OF poc,? 


Quality 
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or.c/.:’.'-- - _ v 

OF POOR Q‘- ‘ * 


Standard triangle 


The BSG idea: Loft & Boolean sum 

P 1 F " Vl^q^ 0 '^ + ^(OFCl-q.q) 

+ ^(OCl-q^ Q (0,q) + K 1 (O(l-q)F 10 (l-q,q) 


( 0 , 1 ) 




use modem version in 
practice. 



a,8,Y polynomials from the Birkhoff Pit. 

Convex Combination otT.,F + ST2 F + YI3F 
T^F C^ interpolant cn edge i 

A A A 

Brown a.S.Y rational functions from Shepard's Formula. 

Little 




Rational, almost polynomial. 
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OnlS*N«- 

OF POOR 


q'JAUTY 


(2) Clough-Tocher 


piecewise cubic 




o 

Problem: Find C Clough Tocher. 

Tool: Farin' s Bdzier Triangle methods 1979. 


Solution: Barnhill, Farin, and Little 1980. 
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Contouring: This is sometimes the problems, e.g., hidden surfaces, silhouette 

edges . 


Adaptive subdivision schemes 


Little 

Petersen 
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Point Methods : Arbitrarily Spaced Data 


Shepard's Formula (1968) 


(SF)(x,y) 


y F i 

l ~7L 

1 d i 

1 -i- 

i d t 


(x,y) 7 s (x^y^ 


(x,y) = (Xj,y) 


Solar maps 


where d^ « d^(x,y) « distance from (x,y) to (x^ , y^) 

l V!k ^ 

Rewrite SF « ■ » = T w.F. 

TTiTS i 11 


Wi 


where v^Xj .y^) 


1 if i - j 

0 if i j* j 


cardinal form. 


SF interpolates and is continuous. 
Global method / flat spots 


Inprovemcnts : Barnhill & Poeppelmeier 1975 
Franke 1975 
Vittitcw 1978 


Little 1978 
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Remarks : Patch methods are local methods. 

Shape fidelity requires at least local quadratic precision. 

Interactive design - real time computations. 

Interactive viewing - use the hardware. 


Reference: R. E. Barnhill, Representation and Approximation of Surfaces, 

Math. Software III, J. R. Rice, ed. , Academic Press, 1977. 
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variables . 

Tessellation of 3D domains into tetrahedra. Little 

I 

4D Surface Interpolants Gregory 

| Mansfield 


3D Contours 
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BACKGROUND 

This paper is devoted mainly to a physical and geometric interpretation of 
the surface fitting technique discovered by the author in 1963, which was 
called multiquadric equations in 1971 (Hardy 1971). It was not until 19S0, 
after a report by Franke (1979) that it was recognized that multiquadric 
equations or MQ may be interpreted very simply as a linear combination of 
three dimensional distance functions. The similarity of the HQ method to 
a simple summation associated with point mass models in qeodesy was rec- 
ognized very early (Hardy 1972). As it turns out now this was simply the 
other side of the biharmonic-harmonic "coin", ard I have recently coined 
a multiquadric label for the point mass anomaly also, i.e. reciprocal 
multiquadric or MQ . In this case we construct a set of point mass 
anomalies having positive and neqative values as contrasted with always 
positive masses. By requiring or assuming the sum of mass anomalies to 
total zero we do not change the total mass and are therefore dealing with 
irregularities in a distribution of mass with respect to what we perceive 
to be some standard distribution of mass. For disturbing potential outside 
the anoialous masses we obtain a solution with a linear combination of 
three dimensional reciprocal distance functions, in which the originally 
unknown point mass anomalies are treated as undetermined coefficients. A 
reciprocal distance function is harmonic, satisfying Laplace's differen- 
bal equation. Hence, a linear combination of such functions is also har- 
monic. An alternative way of looking at the problem and its solution is 
to consider the integral 

T(r p*W = fi'V/f-' l (r .0 ,X :r,fl,X) dm 

in whicn p is a point on or outside the spherical body where disturbing 
potential T is measured. This intearal cannot be formally integrated 
because the equality dm = fi(r,0,X) dv contains an unknown density func- 
tion 6(r,i*,X) under the integral sign. Therefore the integral, considered 
in this light, has the characteristics of an integral equation. Jaswon 
and Synm (1977) have studied problems of this type in both potential and 
elasticity. It is this form of a numerical approximation to a linear 
inteqral equation which provides a solution for the anomalous density 
function using ’X)' . fieasurements of disturbing potential at n points 
and the formation of a system of up to n linear equations provides the 
f. ndanental basis of this approach. However, the approximation of the 
density function is not the onmary goal usually. After obtaining a good 
approximation of the density function it is used in the su'imation form 
to evaluate T at any point, usually where F has not been measured. The 
ultimate outcore then, is a procedure that one can classify as prediction, 
approximation, or surface fitting, if not some other form of numerical 
analysis. 


A frequent problem with HQ* , or point uass models in general, is to rind 
the best depth or ladius for placing the anomalies. Ha*'dy (197C, 1979) 
and Hardy and Gopfert ( 1976) have provided a very sati 1 factory solution for 
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this with respect to both spheres and planes. This is the "best-r formula”, 
but it does not seem to be well known. The original MQ form (now known to 
be biharmomc) is relatively insensitive to this problem, and as noted by 

Franke (1979) it seems to give results equal to or better than MQ’ 1 for 
most surface fitting purposes. Hence HQ essentially as used in 1971, will 
probably become the favorite of the tv/o surface fitting techniques. 

The emphasis below will be on MQ rather than MQ" 1 but contrasts and similar- 
ities of the tv/o will be shown. As will be seen, it, is the geometric and 
physical interpretation that has been applied to MQ -1 which has contributed 
to a better understanding of MQ. To summarize briefly in advance: 

(1) MQ’ 1 is harmonic; MQ is bihamomc. 

(2) MQ’ 1 deals basically with exterior disturbing potential and 
satisfies Laplace's equation; MQ deals basically with interior 
and surface displacements, elastically, and satisfies Poisson’s 
equation. 

(3) In both cases the solutions may be viewed as beinq numerical 
approximations of an integral equation in which an unknown den- 
si ty. function is the physical source for disturbing potential 
(MQ L ) or elastic displacement (MQ). 

Most of what follows lias been taken quite literally from my rect-H papers 
(Hardy 1980,1981). There has not been a rapid expansion of my knowledge 
on the subject since 1980; however, ! am taking advantage of this opportun- 
ity to renedy a few misleading statements, to correct outright mistakes, 
and to chanqe other matters, particularly Figure 2 which represents the 
elastic displacement of a sphere. Hence what is presented is considereo 
to be a modest improvement over my previous papers. 


B I HARMONIC-HARMONIC MODELS 
FOR SURFACE FITTING 

Recognition that MQ is biharmomc in three dimensions, just as MQ’ 1 is 
harmonic in three dimensions, was expedited by Franke’s (1979) descrip- 
tion of Duchon's thin plate spline or TPS. Franke noted similarities 
of IPS and MQ in the fact that ordinates for a single kernel function 
get larger in both cases as the distance increases. TPS involves a 
biharmomc function in two dimensions of the form: 

r log r with r = (x”i-y ) : 

whereas MQ involves a biharmomc function in three dimensions of the form: 

2 -1 2 2 2 
r *r or simply r with r - (x +y +z ) *. 
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Host references in the theory of elasticity deal more with two dimensional 
theory than that of three dimensions. Nevertheless Inharmonic functions 
of both types show up in mathematical physics, particularly in cases in- 
volving relationships of ootential, elasticity and hydromechanics. How I 
will show you several groups of equations and give brief comments on each 
group. 


n 


MQ MODEL: 

Z da Q. = H 

J=l J J 

0) 


q j = [ (x - v 2 + (Y ■ v 2 + s2 ] 1/2 

(2) 

DATA EQUATIONS: 

n 

Z da Q • = H. l = 1,2,. ...n 

J=1 - 1 

(3) 


1 




* (V,- y 



(4) 


In this group we see the original MQ method. Ordinates of 
H consists of linear combinations of hyperboloids centered 
at data points. 6 was considered as a constant, whereas 
we will see later that it can be treated as the difference 
between a constant Z-| and a variable Z 


MQ" 1 MODEL: 


G Z da Q" 1 = T 
J=1 J J 


Q ” 1 = f (x - X ) 2 + {Y - Y ) 2 
J J J 



(5) 

( 6 ) 


DATA EQUATIONS: 


n^ 

G ) da 
1 = 1 J 




i = 1 ,2,. . . ,n 


(7) 


' [ ( \ ' X / * < V , - \ 1 ? * < 2 , 



ibe reciprocal MQ model above is actually a point mass 
anomaly model for disturbing potential T Q to the minus 1 
is a continuous reciprocal distance function m three vari- 
ables. For computational convenience we can locate point 
mass anomalies at a constant depth 5 = Z, . We can also 
make all measurements of T on the XY plane at Z =* 0. Then 
(Z-Zj) becomes the 6 in the MQ equations of the previous 
group 



r _ r -1 T r 1 ORIGINAL PACE 13 

SOLUTIONS: [daj =■ [Q^J" [Hj OF POOR QUALITY 

h] - M 00 

PREDICTIONS: [qj [qJ 1 [«,] - [» p ] (U 

•Mkjr’w-w (u 

The solutions and predictions for MQ and reciprocal KQ, as 
above, follow the same basic pattern in each case. We 
don't need to consider these details now since this is not 
the main purpose of this paper. 

LET 5 2 - (Z - Zj) 2 ' < 13 


Q. = [(X - x/ * (V - x/ ♦ (Z - z/]’ 


The equivalence previously mentioned and identified above, 
causes the MQ basis function to be a Cartesian distance 
function in three variables, analogous to the reciprocal 
distance in three variables . The implication is present 
then, that the undetermined coefficients daj in the two 
cases should have the same physical moaning. This is veri- 
fied by the mathematical theory of elasticity. 

KQ MODEL ( 81 HARMONIC) 


l 2 (7 2 Q) « ? 4 Q = M v + 


4 

+ 


3Y t 3Z‘‘ 



(16) 


Any single 0 or distance function as above satisfies the 
biharmonic differential equation m three variables as 
given. Thus a linear combination of all distance functions 
used in MQ approximation is biharmonic. 
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J 2 fQ - 1 

3Z 2 




07 ) 


(18) 


q-1 as above, is the generating function for zonal harmon- 
ics, and through die decomposition formula, leads to com- 
plete spherical harmonics satisfying the Laplace differen- 
tial equation. Thus a linear combination of reciprocal MQ 
functions satisfies Laplace's equation. 


DUCHON'S TPS: 


w(p) = I? - qi 2 log Ip - q! 


(19) 


!i> - ii - [(Xp - x,) 2 * t’p - '/] 


1/2 


( 20 ) 


The basic idea of Duchon's TPS is given in equations (19) 
and (20) above v.’(p) is the deflection at tne point loca- 
tion (XnYp) where a concentrated load A is applied. (Xq, 
Yq) is a point located on the boundary defined as a simple 
support for the plate. D is the constant of structural 
rigidity. Then Ip-ql is the distance between 2 points in 
the same plane. I'ote that if we let Ip-ql = r then equation 
(19) is in the simplified form kr- log r. Log r is the 
well kno.vn loginthmic potential m 2 variables. 

D’UCnCVS MODEL * 





(Y - Y, 



♦ a, A + a - N + a- 1 f(X,Y) 

I <- J 


(?D 


Duchor's TPS model given above indicates that f(X,Y) is a 
linear combination of terms r2 log r plus three terms that 
physically account for rigid bouy displacements. These do 
not affect stress or strain m a thin plate 



— i tl"* * r 




1 
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£ A J f (X < ' */ * (V i " Y /] 103 [ (X ' ■ V 2 * (Y < ' V 2 ]' /2 

+ a 1 X 1 + a 2 Y i + a 3 = f (X 1 ,Y^ ) (1 = 1,2 n) 


WITH tONOITI ONS : 
n n 


£ A ’ 0; £ A.X . « 0; £ A.Y » 0. 
j-1 3 J-1 J J J=1 J J 


( 22 ) 


The above equations show that n measurements of stress or 
displacement ordinates f(Xi,YiT are sufficient to determine 
n valuer- of upper case A coefficients (concentrated loads), 
when 3 condition equations involving rigid body motion are 
included in the '.ystem of equations. 

ONE FORM OF HARMON I C-B l HARMONIC RELATIONS 

X = >' 2 $> + Y (23) 

X , B I HARMONIC; i AND?, HARMON I C . 
r, DISTANCE IN 2 OR 3 VARIABLES 


This relationship given above is one of several classical 
statements concerning biharmonic-harmcnic relations in the 
theory of elasticity. It is applicable to both HQ and TPS. 

CJCHOrS TPS BASIS 

1/2 

X = r Z LOG r + Y r = (X 2 * Y 2 ) " (24) 

V 4 X * 0 BIHARMWIC IN TWO VARIABLES 


Duchon's TPS is bihamonic m two variables because it can 
be expressed as a combination of two namonic functions, 
one of which is log r, harmonic in the two variables , multi- 
plied by the square of r in two variables. The other 
function is C- “ a X + a Y + a , also harmonic in two vari- 
ables. 12 3 

HARDY'S MQ BASIS 


> = r 2 (1/r) ♦ Y » = ( v - 2 i ) Z * 2 2 ) 


1/2 


V X 51 c 


£ : HARMON I C IN T^REt VARIABLES 


( 25 ) 
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MQ is biharmonic in three variables because it can be ex- 
pressed as a combination of two Harmonic functions, one of 
which is MQ"* m three variables, multiplied by the square 
of r in three variables. The other harmonic function will 
be discussed later. 


Figure 1 illustrates the cross .section of a sphere after 
transformation of MQ and reciprocal I IQ to spherical coordin- 
ates. We assume the sphere is an idealized solid elastic 
body with constant density or density as a function of the 
radius only, except for one mass element. A mass excess 
of daj at a single element induces a disturbing potential 
of the otherwise sohencal cquipotent lal with respect to 
the sphere. The magnitude is greatly exaggerated to show 
clearly the shape of the disturbed equipotential surface. 
MQ - 1 approximation as a whole consists of a linear combina- 
tion of ordinates of such disturbed surfaces. 


original r\;~ 
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FIG 1 DISTURBING POTENTIAL 



r 
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FIG. 2 BIHARMONIC DISPLACEMENT 


Figure 2 is the same diagram except that it shows m exag- 
gerated form the biharmonic displacement of the solid spher- 
ical body caused by the same point mass anomaly d a j . MQ 
approximation as a whole consists of a linear conDinaticn 
of ordinates of such disrurbed elastic surfaces. Note that 
the displaced surface is inside the spherical surface for 
a positive anomalv, due to physical contraction of the 
body. The neqative displacement is least nearest the anomaly and increases 
negatively as the spherical distance increases. An important point to be 
visualized with this illustration is that the point mass anomaly ch is 
itself displaced radially inward (for a positive anomaly) durma the inter- 
action of the standard masses and the point mass anomaly. This change in 
position induces a small chanqe in the exterior disturbino potential. 

This induced chanqe in the exterior harmonic function is probably associated 
with the form of the harrr,onic-bi harmonic relations illustrated earlier. 

In brief there appears to be justification for addino three terms in one 
variable each becajse of three dimensional displacements in the solio body, 
and possibly a constant term as well to completely fulfill the theory of 
elasticity. Practically it doesn’t seem necessary to include these terns 
when the MQ nethod is aoplied to non-elastic problems, as a general approx- 
imation scheme. I suspect the statement would be true for buchon’s TPS. 



CONCLUDING REMARKS 
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I wish to comment briefly on a possible reason why the MQ 
methods gave generally better results than Duchon's TPS m 
Franke's study. Duchon's method involves direct applica- 
tion of externally concentrated forces at the surface of a 
reference plane; there are no body forces. The MQ methods 
use body forces induced by anomalous gravitation; there are 
no concentrated external forces. Hence the MQ biharmonic 
function is generally a smoother function than Duchon's TPS. 
This may account for some differences in the approximation 
properties of the two methods. 
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BSPLASH: A Three-Stage Surface 

Interpol ant to Scattered Data 


By Thomas A. foley 

California Polytechnic State University 
Computer Science Department 
San Luis Obispo 


ABSTRACT 

Given N distinct points ( x - » y.) and N real numbers z^, 

BSPLASH constructs a function G (x, y) that satisfies G (x^ y i ) - z., 
tor i = 1,..., N. This C 2 ir.terpolant consists of a bicubic spline 
approximation and Shepard's bivariate mterpolant. 
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1. Introduction 

This paper presents a three-stage procedure that solves the 
following bivariate interpolation problem. Given N distinct points 
in the plane(x.j, y^) and N real numbers z^, construct a function 

G(x,y) that satisfies G (x.,y.) = z* for i * 1 N. This is re- 

ferred to as the scattered data interpolation problem because the 
data points (x,. ,y.) are not assumed to fall on a rectangular grid. 

The interpolation problem can be interpreted as fitting a sur- 
face through N points in three dimensional space and thus has many 
applications. In mineral exploration, exploratory wells are drilled 
and the depths of various layers are recorded. Given this data, 
surfaces representing these layers can be constructed using inter- 
polation methods. Such an example was studied by Robinson, Charles- 
worth, and Ellis [9] in petroleum exploration. Foley [3] and [4] 
used bivariate interpolation in the characterization of radio - 
nuclide activity resulting from nuclear tests. Samples of activity 
were measured at various locations, the (x.,y_. ) points, and the 
magnitude of the readings were the z.. ' s . The survey paper by 
Schumaker [10] gives applications in medicine, computer aided de- 
sign, electronics and geology. 

The new approach presented here is similar to the delta 
iteration methods of Foley [3] and Foley and Nielson [5], but the 
new method is more stable, visually smoother on smooth data, and 
it uses less storage. It is equally efficient on large data sets 


jfstt 
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and is capable of smoothly filling in arpas that are void of data. 
This globally defined intcrpolant can be displayed by a three- 
dimensional surface, a contour map, or by a table of interpolated 
2 values. 

The name of the algorithm is called BSPLASH because it uses a 
bicubic spline approximation and a modified Shepard's intcrpolant. 
The motivation for this approach is based upon the fact that many 
methods that apply directly to scattered data do not yield smooth 
or desirable surfaces. On the other hand, many methods that yield 
efficient smooth interpolants only apply to data that fall on a 
rectangular grid. 

2. Modified Shepard's Method 

It will be convenient to use operator notation and to assume 
that there is some underlying function f(x,y) defined at the data 
points that satisfies ffa^.y^} - z- for i = 1,...,N. An easily 
implemented scattered data interpolant is the rodified Shepard's 
method described m Foley [3] that is defined by 
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r 


N 

E 

S[f] (x,y) = < i=l 


N 

v* 

t-i 

1=1 


Pi ( x »y) 


for (x,y) f (Xj.yj) 
j — 1 , . . . , N 


Pi (x,yT 




for (x.y) = (x^.yj). 


where d i = (x - x.) 2 + (y - y.) 2 , r. is the distance squared from 
(x., y^) to its 5th nearest data point divided by four, and 

d i< R i 4 


D.. (x.v) = 



'Wi-rv ’ r ?" ; ? -rr.V¥ 
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This method produces an in^erpolant based on the inverse of the 
distance from a point to the data points. The proof of the following 
theorem can be found in Foley [f] and in Gordon and Wixom rj3], 

*■ 7 

Theorem 1 

Given N distinct points (x^.y^), then 

a) S is a linear operator; 

b) S[f] (x.,y.) = f (x^) for i = 1 N; 

c) SCO e C°° (R* 1 ), that is S[f] has continuous partial derivatives 
of all orders for all (x,y); 

d) S[f] (x. ,y^) = 0 and |7 S[f] (x^y.) = 0 for i = 1.....N; 

e) S[f] satisfies the max-min principle min f(x.,>.) < S[f](x,y)< 

i<N 

max f(x.,y.) for all (x,y); and ~ 

i<N 1 1 

f) 5[f] is invariant under translations, rotations, and magnifi- 
cations of the data points (x^,y.). 

Properties a), b), and c) state that S[f] solves the scattered 
data problem with a continuous function, while d) says that S[f] 
is flat at the data points. Property e) is important when dealing 
with data that has a large variation in in small regions. S[f] 
will not oscillate violently as some other rethods might. The 
final property implies that S[f] depends on the relative distances 
between data points, and not on the placement of the axes, nor on 
whether distances are measured in inches or miles. 

Modified Shepard's int^rpoiant is very fast computationally, 
requires very little storage, and easily generalizes to higher 
dimensions. Unfortunately, even though S[f] c C™ the plots are 
not visually smooth. 

CS:c:r*7L 
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Figure 1 displays the plot of the bivariate function 
(1) f(x,y) = .75 cxp ( - (9x-2) 2 * _(9 y-2l 2 } + - 75 exp(-l9x+i) 2 - 9y+l) 

+ -5 exo(- (9x-7) 2 + (9y-3) 2 l - .2 exp(-(9x-4) 2 - (9y-7) 2 ) , 

4 

on the domain 0 £ x _< 1 and Ojc y £ 1. 

From this surface, 100 data points were chosen whose (x.y) 
coordinates are shown in Figure 2. Each point was chosen randomly 
from a uniform distribution on a square with side length 1/9 centered 

at (i/9, j/S) , i,j * 0,1,2 9. This function and data are used 

here because Franke [6] used them as his primary test case in com- 
paring many bivariate intcrpolants. 

Figure 3 shows the modified Shepard's interpolant applied to 
this data. The maximum absolute error is .2403 and the average 
absolute error is .0284. These enors were computed using the 
differences at the 33 by 33 grid used to plot the surfaces. The 
values range from .027 to 1.17. 

3. Bicubic Splines 


Another major corrporent of BSPLACH is the bicubic spline that 
solves the following gndded data interpolation p.ohlem. Given 
(XG , YG , ZG. .), i - 1,...NXG and j - 1....NYG, construct a tunction 
H(x,y) that satisfies h(XG i ,YGj! = 76^. The corresponding operator 
notations assumes that there is some underlying function g(x,y) 
defined at the grid points that satisfies gJXG^.YG^) = aG^. 
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The points (XG. ,YG ) will be referred to os the rectangular grid 
• J 

points. 

The natural bicubic spline interpolant to g(x,y) over the grid 
(XG..YG.) is denoted by D[g] (x,y) and it solves the gridded inter- 

• J 

polation problem B[g] (XG {> YG>) = g(XG-.YG-) for i 3 1,...,NXG and 
j 3 1 . ,NYG. 

B[g] is formed by the tensor product of natural cubic splines 
in the x and y directions. The function is a piecewise bicubic poly- 
nomial of the form 

3 

E E b u * 

1=0 j=0 

on each rectangle determined by the rectangular grid points (XG-.YGj). 
It is pieced together smoothly so that B[g] has all of its second 
order partial derivatives continuous. B[g] also minimizes the curva- 
ture functional : 

L(f) -J ft ?i 2^ x »y)) z dx d y 

over all functions f(x,y) that solve the same gridded interpolation 
problem and satisfies certain continuity conditions. See deBoor [1] 
for a detailed description. 

Other bicubic spline inter polants exist that use different end 
conditions. The natural end conditions were used here primarily be- 
cause they were easily accessible in the software package I.M.S.L.[8] 
in the subroutine IBCIEU. This subroutine was »?plied to the test 
function (1) using an equally spaced 9 by 9 grid on the unit square 
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and this is shown in Figure 3. The maximum absolute error is .0406 
and the average absolute error is .0033. 

In many cases, this interpolant is a good choice for gridded 
data because it is computationally efficient, visually smooth, it 
belongs to C , it is globally defined, and it is easily accessible. 
Unfortunately, it only applies to gridded data. 

4. BSPLASH 

BSPLASH uses a bicubic spline and modified Shepard's method 
in the second and third stages of the procedure. The first stage 
of the alognthm is to generate a gridoed data problem. This con- 
sists of defining the rectangular grid points (XG.,YG-) and the 

* J 

values ZG. = g(XG ,YG ) using local least squares approximations. 

• J » J 

The second stage is to form the bicubic spline through these points. 
The final stage adds a correction term using Shepard's method to the 
bicubic spline so that the scattered data interpolation problem is 
solved. 

BSPLASH allows the user to enter his own rectangular gridpoints 
or else it ccrputes the grid for him. The grid algorithm is first 
applied to the x 1 ,s and then to the y. ' s . The objective is to 
compute grid points that cover the data points (x^.y^) proportion- 
ally to the density of the data points without having grid points 
too close together or too far apart, and to have all the data 
points fall inside the grio. There is a restriction that NXG 25 
and NYG £ 25 for large data sets. 

Let H = IR0UND(/N), k = IR0UND(N/M) , NXG = K+2 and NYG = K+2. 
Sort t w • x-coordinates into increasing order. Set XC 2 to the average 
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of the first k x-coordinati ., XG-j to the average of the next k small- 
est x-coordinates,..., and XG^ to the average of the k largest 
x-coordinates. 

Let U = {XG^ ■ XG2)/(H-1). If the H grid points were equally 
spaced, then U would be the difference between two consecutive grid 
points. Set XGj = Xj - U and XG^ = x^, + U so that all of the data 
points fall inside the grid. 

While the interior grid points are being computed, consecutive 
grid points are compared to see if their difference is between U/2 
and 3*U. If their difference is less than U/2, they are considered 
to be too close and they are averaged together thus reducing NXG by 
one. If their difference is greater than 3*U, a new grid point is 
inserted at their midpoint and NXG is increased by one. 

The y-coordinates of the rectangular grid are defined in the 
same manner. 

Figure 5 shows the results of applying this grid algorithm to 
sets of data consisting of N - 100, 33, and 25 points. The grid 
points are where the orthogonal lines intersect. Note that the 
second highest horizontal line in Figure 5b is the average of two 
grid lines that were too close together. All three of these (x,y) 
data sets were used by Franke [6] in his comparison of many bi- 
variate interpolants. 

The rest of the first stage is to define the values ZG^ = 
g(XG ,YG ) at the grid points. For each of the grid points (XG , 

' J I 

YGj), find the seven nearest data points (x^,y^). To simplify the 

notation, assume that the nearest data points to (XG ,YG ) are 

^ J 
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(x i ,y l ),...»(x 7 ,y 7 ). Let d fc = (x^-XGj) 2 + (ytf YG j) 2 « and note that 
dj 1 ^2 — " - ^ 7 * A weighted least squares fit to (x^.y^.z^,), 
k = 1.....7, by a quadratic is formed and then it is evaluated at 
(XG-.YG ). That is, solve 

* w 

(2) rain £ J- ( 9 <x k ,y k ! - z/ 

a. k=l K 

J 

2 o 

for the Sj ' s , where g(x,y) - + a 2 x + a^y + a^x + a^xy + a^y . 

Then define 7.G.. = g(XG. ,YG-). 

' J ' J 

To add stability to this process, BSPLASH sets ZMItJ = min 
(z, ,. . . ,z 7 ), ZKAX = max (z^,...,z 7 ), and then defines the grid values 
by 

{ ZMIN if g (XG^YGj) ^ ZfUK 

ZKAX if g (XG,,Y6) > ZKAX 

* J 

gJXG^.YGj) otherwise 

The I.H.S.L. subroutine LLSQF is used to compute a^,a 2 ,...ag 
for each of the grid points. This subroutine will properly handle 
the case where the minimization problem (2) has many solutions by 
using the fit of lowest degree. 

This first stage can be used to define function values at any 
point, but since this depends on the seven nearest data points, the 
function may not be continuous. However, this weighted least 
squares approach gives a good local approximation to f(x,y) at the 
rectangular grid points. 
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The second stage is to form the bicubic spline B[g], where 

g(XG., YGj) = Z G.y 1 « 1 NXG, and j * l,...,fiVG. Although this 

generally yields a smooth approximation to f(x,y), it does not inter- 
polate the scattered data (x-.y^.z.). 

The third and final stage solves the scattered data problem by 
adding the correction term S[f - B[g]] to B[g] to yield the BSPLASH 
interpolant. 

P[f] (*,y) =- S[f - 8[g]] (it ,y ) + B[g] (x,y). The correction 
term uses the modified Shepard's method to interpolate the differ- 
ences between z. and the bicubic spline B[g] evaluated at (x ,y ), 
i = 1 . ,N. 

Theorem 2 

Given N distinct points (x^.y^) and N real numbers z^ = ffx^.y ), 

then 

a) P[f] solves the scattered data problem P[f] (x ,y.j) = fjx-.y^) 
for l = 1 , . . . ,N; and 

b) P[f] has all of its second order partial deri atives contin- 
uous for all (x,y). 

Proof: 

By b) of Theorem 1, it follows that P[f] (x-.y^ 

= STf - B[g]] (x i ,y 1 ) B[g] (x^.y^ 

= f^.y,) - ^[g] (x^y,) + B[g] ( x ^ ,y i ) 

B f (x i ,y i ) . 

Since P[f] is the sum of C*" Shepard's correction function and the 

? ? 

C bicubic spline, it follows that P[f] belongs tc C for all (x,y). 

Q.E.D. 
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Figure 6 shows the results of BSPLASH applied to the data sets 
of N = 100, 33, and 25 points from (1) evaluated at the (x,y) points 
shown in Figure 5. The maximum absolute errors are .0443, .2293, 
and .1220 respectively, and the average absolute errors are .0060, 
.0435, and .0277 respectively. 

These discrete errors compare favorably with the best of the 
methods tested by Franke [6]. The visual smoothness of BSPLASH would 
also rank high with those methods. BSPLASH was applied to three 
other Lest function of Franke [6] on the same three data sets of 
100, 33, and 25 points, and the results were accurate and visually 
smooth. 

The storage required for BSPLASH is very low. Besides the 
Storage of the data points (x.. ,y^ ,z . ) , i * 1,...,N, less than 3N 
locations are needed to store the grid, the grid's z-values, and the 
local parameters r u$ec in modified Shepard's method. Host tri- 
angular based interpolants reouire storage on the order of 30N and 

2 

some others require storage of more than N elements. 

The execution times can't be accurately compared because dif- 
ferent computers were used. The results here were done on a Cyber 

2 

CDC 170/730. The overall computation tine is on the order of N , 

but the observed times appear linear in N even when the grid algo- 

oM 

rithm was used. The executing:/ times for BSPLASH when N = 25,50,100, 
200, 400, and 800 points wci<> used were 4.5, 8.3, 17.2, 39.6, 104.5, 
and 206.2 seconds respectively. Some other bivariate interpo 1 3nts 
are on the order of f r and are not efficient for large M. 
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APl’.AiDIX 

IlSPLAL'iI on Drankc's Test Data 

DSI'LASI! was applied to several sets of data that I'rankeCG] used in 
the co-.iparison of rtany bivariate interpolonts . Section 1 describes the 
data sets, octsiiBn d - -3 ! -< / *<•*»■ - ■■ frh- e - <1 i >;(.- v i> T'«Trm. cor s 

&&&&, section 3 contains the plots of the 3-D surfaces, and section 4 
lists the discrete errors of the better ref ho Is teste , by i'ranhe. The 
first and fourth sections are edited xeroxes of franco's technical 


report 
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1.2.2. The Test Problems 

The basic set of test problems consisted of six different test functions 
over three different x-y point sets, and two x-y-z point sets from the 
literature, one of those used in a second version with one of the coordinates 
scaled. Another interesting test was the computation of a "cardinal" 
function obtained by setting all function va'iues on a point set to zero, save 
one. 

The six test functions were all to be approximated on [0, 1] . Four of 

them were basically obtained from McLain's paper [39], but were translated 
2 2 

to [0, 1] from [1, 10] and some modified slightly to enhance the visual 
aspects of the surface. The other two were generated by the author to provide 
a fundamentally different shape in one case (saddle), and to provide a surface 
with a variety of behavior on one surface to serve as a principal test func- 
tion. 


The principal test function is given by 

l 2 


f,(x, y) - .75 oxp[ - lasaifistaf] * .75 „p£ - - 2ft!] 


• + 


.5 exp[ - ] _ .2 exp[ - (Sx-4) 2 - (9y-7) 2 ]. 


This surface consists of two Gaussian peaks and a sharper Gaussian dip 
superimposed on a surface sloping toward the first quadrant. The latter 
was included mainly to enhance the visual aspects of the surface, which is 
Shown in figure 4. 0.1.0. 

The second test function, essentially obtained from McLain is 
f 2 (x. y) = |{tanh(9y - 9x) + 1]. 

2 

This surface consists of two nearly flat regions of height 0 and g, joined 
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by a sharp rise, alirost a cliff, running diagonally from (0, 0) to (1, 1). 
The test surface is shown in Figure 4. 0.2.0. 

The third test function was generated by the investigator and is 


f 3 (x, y) = 


1 .25 + cos(5.4y) 
5[1 + (3x - l) 2 ] ' 


This surface is saddle shaped and is shown in Figure 4. 0.3.0. 

The fourth test function, essentially obtained from McLain, is 


F 4 (x, y) = ~ exp[ - y^( (x— jj-) 


(y-p )]. 


This surface is a Gaussian hill which slopes off in rather gentle fashion in 
2 

[0, 1] . It can be seen in Figure 4. 0.4.0. 

The fifth test function was also essentially obtained from McLain and is 

2 2 

f 5 (x. y) = lexp[ - §l((x-I) + (y-1) )]. 


This surface is a steep Gaussian hill which becomes almost ; C ro at the bound 
aries of the unit square. It can be seen in Figure 4. 0.5.0. 



The first set consisted of 1 CO points generated by a pseudorandom nunber 
generator, one point m each square of side g- centered at (-g-, g-) for 
i, j = 1, .... 10. This yields a set of scattered points forced to have 
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somewhat uni form density, although as can be seen in Figure 0.1. 0.0. there 
are locally large variations in density. The triangulated set of points is 
also shown in Figure 0.1. 0.0. Part of the unit square is outside of the 
convex hull. The points arc listed in Table 1. 

The second set of data consists of 33 points and was generated by tne 
investigator to purposely have some areas sparsely populated by points 
while other areas are not. Tins set of points is shown in Figure 0.2. 0.0. 

The points are listed in Table 2. J 

The third set of points was digitized by Gregory M. Nielson and is 
similar in disposition to a set of points appearing in McLain [40]. This set 
of points is shown in Figure 0.3. 0.0. Part of the unit square is outside the 
convex hull. The points are listed in Table 3. 

Two sets of data were obtained t'ron the literature, and one of these was 
scaled in one variable to obtain another. A fourth set was used to generate 
a ''Cardinal Function". The da*a given in Table 3, and shown in Figure 0.3. 0.0. 
was given the following function values: ffx^,. /^) - 0 except : 

f(.1875, .2625) = .2. Here .2 was used for visual purposes rather than 1 as 
would ordinarily be done for a true cardinal function. This gives some infor- 
mation about the influence of one point cn the surface for moderate sized 
point sets. Of the two sets of points from the 'literature, one is from Akima 
[1] and was obtained during a study of waveform distortion. It is repeated 
here in Table 5, and shewn in Figure 0.5. 0.0. The second was obtained from 
Ferguson [14] and is repeated here in Table 6, and shown in Figure 0.6. 0.0. 

The same set of data, but with the y coordinate multiplied by three was 
also used to show effects of scaling only one variable, and is shown in 
Figure 0. 7.0.0. For visual purposes, the function values given in Table 2 
are actually .6 more than given by Ferguson. As can be seen from Figure 
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Fl ^oopt-s 


j: Franke - 3 

4: AHra 
jq ; AMiU Mod. I 
1 3 ; Nielson - franke .) 
14: Mod. Quad. Shepard 
15 ; Akira Mod. Ill 
2i: Franke - TPS 

28: Lawson 

19: Nielson MinNorm 

21: Hardy Quadric 
?3: Duchon TPS 
27: Hardy Pecip. Quad. 
30: Foley III 
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Maxi nun 

Mean 

RMS 

Deviation 

Deviation 

Deviation 

.0919 

.00842 

.0148 

.0647 

.00787 

.0125 

.0856 

.00784 

.0133 

.0782 

.00741 

.0122 

.0573 

.00785 

.0128 

.0520 

.00729 

.0117 

.0940 

.00387 

.0164 

.0951 

.00783 

.0124 

.0492 

.00537 

.00940 

.0225 

.00181 

.00357 

.0518 

.00525 

.00547 

.0247 

. f >0233 

.00518 

.0636 

.00473 

•C0941 

.o*{4 3 . 

*0066 


Maximum 

Mean 

RMS 

Deviation 

Deviation 

Deviation 


franke - 3 
AMra 

Aki r.a Mod I 
forlson - Frarke Q 
Mod Quad. Shepard 
Akira Mod. Ill 
Franke - TPS 
Lawson 

Nielson Min’lorn 
Hardy Quadric 
Duchon - Tp S 
Hard/ Pecip. Quad. 
Foley HI 
£ SPLASH 


Method "V", 

1: Frante - 3 

4: Akira 

10: Akira Mod. I 

13 : Nielson - Franke Q 

14; ‘■lad. Quod Shepard 
16. Akira Mod III 
24; Franke - TPS 
28: Lawson 

19. Nielson M.in'.'ora 
21: Hardy Quadric 

23 Duct. on TPS 
27: Hardy Pecip. Q^ad. 

30: Foley III 


1: 

4: 

10- 

13 : 

U. 

16 : 

24: 

23. 

19 : 

21 : 

23 : 

27: 

30 : 


.347 

.158 

.197 

.150 

.184 

.164 

.218 

.287 

.150 

.137 

.153 

.140 

.296 


Maxi rum 
Deviation 

.240 

.134 

.129 

.153 

.153 

.155 

.129 

.202 

'.124 

.119 

.121 

.119 

.165 


.0477 

.0384 

.0400 

.0326 

.0340 

.0372 

.0346 

.0462 

.0305 

.0181 

.0233 

.0153 

.0350 

,o4-2>£~ 


Mean 

Deviation 

.0359 

.0282 

.0280 

.0350 

.0353 

.0355 

.0267 

.0327 

.0235 

.0235 

.0253 

.0214 

.0195 

, 02-77 


.0732 

.0535 

.0570 

.0455 

.0478 

.0521 

.0517 

.C557 

.0437 

.0269 

.0421 

.0244 

.0546 


RMS 

Deviation 

.0485 

.0335 

.0333 

.0478 

.0435 

.0484 

.0374 

.0458 

.0323 

.0322 

.03**8 

.0294 

.0310 
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1: 

Franke - 3 

.0518 

.002CG 

.00530 

<: 

Akitna 

.0520 

.00303 

.00609 

10: 

Akimc Mod. I 

.0473 

.00257 

.00512 

13: 

Nielson - Franke Q 

.0721 

.00205 

.0O6P3 

14: 

Mod. Quad. Shepard 

.0468 

.00264 

.00551 

16: 

Akina Mod. Ill 

.0953 

.00293 

.00800 

24: 

Franke * TPS 

.0295 

.00243 

. 00183 

28: 

Lawson 

.0280 

.00221 

.00448 

19: 

Nielson HinNorn 

.0424 

.00181 

.00134 

2i 

Ha>-dy Quadric 

.0244 

.00177 

.OCm30 

23. 

Duchon - TPS 

.0344 

.00210 

.00436 

27: 

■ Hardy Recip. Quad. 

.0379 

.00192 

.00338 

30: 

Foley 111 

.0281 

.00223 

.00419 


BSP4-/1SH 


.coal 



Deviations fron Cliff te c t surface, 100 points 
Tabic 0-1.2 


I: 

4- 

••v 

•i- 

14: 

•t- 

.’ 4 : 

r. 

19 : 

11 : 

n. 


1: 

4: 

10: 

13: 

K: 

16: 

24 : 

23 : 

19: 

21 : 

23 : 

27 : 

33: 


franke - 3 
Akira 

Akira Mod. I 
Nielson - Franke Q 
Mod. Quad. Shepard 
Akira Mod. Ill 
franke - TPS 
Lawson 

Nielson 'iin’iorm 
Hardy Quadric 
Ouchon - S 
Hardy Recip. Quad. 
Foley III 


Table D.2.2 


Franke - 3 
Ak ira 

Aki~a Med. 1 
Nielson - Franke Q 
Mod. Quae. Shepard 
Akira Med. ill 
Franke - TPS 
! iwsenr 

Nielson Mr r, Norn 
Hardy Quadric 
Duchcn - TPS 
Hardy Rccip. Quad 
Foley III 


Bsplash 


,0776 

.0124 

.0190 

.0543 

.00350 

.0133 

.0513 

.00747 

.0122 

,0370 

.0137 

.0219 

.0876 

.0121 

.0205 

.0680 

.0106 

.0176 

.0561 

.00913 

.0147 

.0956 

.0126 

.0205 

.0582 

.00800 

.0140 

.0577 

.0129 

.0170 

,052b 

.00777 

.0134 

.0500 

.00353 

.0130 

,0914 

.0165 

.0262 


eC&lCl 


surface, 

33 points 



— j~ 



J2_ 


.1(71 

.0225 

.0103 

.0999 

.0148 

.0257 

.0987 

.0143 

, 0 ' j 2 

.148 

.0166 

.0304 

.163 

.0166 

.0314 

.146 

.0104 

.0305 

.106 

.0143 

.0257 

.132 

.0164 

.0233 

.0942 

.0133 

.0242 

.0995 

.0143 

.0231 

.101 

.0135 

.0235 

.105 

.0139 

.0236 

.0232 

.0105 

.0250 

♦ 0779 

.oio7 
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Franke - 3 
Akira 

Akimu Mod. 1 
Nielson - Franke Q 
Mod. Quad. Shepard 
Akira Mod. Ill 
Franke - TPS 
Lawson 

Nielson HinNem 
Hardy Quadnc 
Duchon - TPS 
Hardy P.ecip. Quad. 
Foley III 

B£P^S14 


MAX, ME AfJ 1,5 


0193 

.00164 

.00294 

0271 

.00224 

.00423 

0254 

.00193 

.00367 

0163 

.00110 

.00206 

0125 

.00112 

.00194 

0142 

.00105 

.00202 

0165 

.00157 

.00273 

0565 

.00149 

.00359 

0195 

.00091 

.00200 

00461 

.00025 

.00052 

,00597 

.00049 

.00092 

,00928 

.00063 

.00135 

,0117 

.00117 

.00196 

,oi9r 

.C OlO 



Deviations from Saddle test surface. 100 paints 

Table D.1.3 



1: 

Franke - 3 

.111 

.0121 

.0224 

4: 

Akir.a 

.0578 

.0110 

.0165 

10: 

Akima Mod. 1 

.0573 

.0104 

.0156 

13: 

Nielson - Frani e v, 

..0679 

.00939 

.0146 

14: 

Mod. Quad. Shepard 

.0724 

.00907 

.0139 

16: 

Akima Mud. Ill 

.0597 

.0104 

.0162 

24: 

Franke - TPS 

.0662 

.0109 

.0175 

28: 

Lawson 

.0535 

.0133 

.0199 

19: 

Nielson NinNorm 

.0571 

.0102 

.0159 

21: 

Ha. dy Quadr.c 

.0262 

.00442 

.00689 

23: 

Duchon - TPS 

.0574 

.00912 

.0140 

27: 

Hardy Pecip. Quad. 

.0505 

.00571 

.00970 

30: 

Foley III 

.0835 

.00333 

.0148 


B>SP2~ASH 

.073& 

.oioS 



Deviations from Saddle test surface, 33 points 

_rr 


1: Franke - 3 

.9533 

.0111 

.0171 

4 : Ak i " a 

.0861 

.0121 

.0202 

10: Akira Mod. I 

.0366 

.0119 

.0203 

13: NtI son - F rar.ke Q 

.0794 

.0115 

.0189 

14. Mod. Qcad. Shepard 

.0759 

.0114 

.0183 

16: Akira Mod. Ill 

.0787 

.0116 

.0139 

24: franke - TPS 

.9714 

.06933 

.0171 

28: Lawson 

.0375 

.0126 

.0205 

19: Nielson Mir. Norn 

.0704 

.0100 

.0172 

21 : Hardy Quadnc 

.0397 

.00570 

.00952 

23: Duchon - TPS 

.0585 

.00310 

.0137 

27: Haidy Recip. Quad. 

.0443 

.00528 

.00955 

30: Foley III 

.0323 

•00S53 

.0165 



.cOoST 


Deviations from Saddle 

test s"rfacc. 

25 points . V- 
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OrtlG'.rtAL FAC- “'-J 

OF POOR QUALITY f\A (\Y . 


MEAfJ 


1: Franke - 3 
4 ; Akina 
10: Akina Hod. I 

13: Nielson - Franke Q 

14: Hod. Quad. Shepard 
16: Akina Mod. Ill 
24: Franke - TPS 

28: Lawson 


.0114 

.00122 

.00109 

.0101 

.00124 

.00177 

.00675 

.00102 

.00113 

.00517 

.00053 

.00003 

.00388 

.00065 

.00009 

* .00330 

.00049 

.00070 

.00560 

.00103 

.00141 

.00899 

.00061 

.00109 


19: 

21: 

23: 

27: 

30: 


Nielson Minflom 

.00303 

.00047 

Hardy Quadric 

.00102 

.00005 

Duchon - TPS 

.00294 

.0001 7 

Hardy Rerip. Quad. 

.00227 

.00034 

Foley III 

.00604 

.00083 

££PLF\^ 

.0077 

♦ CO0& 


Deviations from Gentle test surface, 100 points 


Table D.I.* 


.00069 

.00011 

.00030 

.00060 

.00117 


1: 

Franke - 3 

.0446 

.00603 

.0101 

4; 

Akima 

.0167 

.004S7 

.00623 

10: 

Akina Hod. I 

.0160 

.00442 

.00573 

13: 

Nielson - Franke Q 

.0312 

.00422 

.00637 

14: 

Hod. Quad. Shepard 

.0272 

.00451 

.00679 

16: 

Akina Hod. Ill 

.0204 

.00394 

.00555 

24: 

Franke - TPS 

.0339 

.006S1 

.0107 

28: 

Lawson 

.0269 

.00652 

.0031 5 

19: 

Nielson MinNorm 

.0214 

.00371 

.00563 

21: 

Hardy Quadric 

.03724 

•00 1 ?! 

. COti o * 

23: 

Duchon - TPS 

.0259 

.004 1 5 

.00714 

27: 

Hardy Recip. Quad. 

.0188 

.00206 

.00185 

30: 

Fole> III 

.0349 

.00133 

. OOo/ -* 


^SPLASH 

-ciSiS 

.604-7 



Deviations from Gentle 

test surface. 

33 points ^ 


1: 

Franke - 3 

.0247 

.00491 

.OC651 

4: 

Aki-a 

.0256 

.00541 

.cc:i5 

10: 

Akina Hod. I 

.0240 

.C0541 

.026S1 

1): 

Nielson - Franke Q 

.0340 

.00563 

.00746 

14: 

Mod. Quad. Sheoaro 

.0227 

.C05?« 

.02669 

16: 

Awma Mod. 1 1 1 

.0232 

.00575 

.00760 

74. 

Franke - TPS 

.0245 

.00140 

.02566 

:i: 

Lawson 

.0234 

.0039? 

.005 1 1 

19: 

Nielson HinMorn 

.0161 

.00307 

.00433 

21: 

Hardy Quadric 

.00709 

.00107 

.00158 

71. 

Djchon - TPS 

.0128 

.00265 

.00351 

77: 

Hardy Rerip. Quad. 

.00528 

.00055 

.002;^ 

30: 

Folev’ III 

.0224 

.004 jo 

.00:55 



.Oc^ai 

.co37> 



Deviations frcn Gentle test surface, 25 points 



Table D.3.4 
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ORIGINAL PAUL ? *3 
OF POOH QUALITY 


MAX » 


ME/; tJ 


1: FranLe - 3 

4: Akima 
10: Akiiua Mod. I 
13: Nielson - Franl.e 0 

14: I'.od. Quad. Sheparo 
1C: Akin-a Hod. Ill 
24 : Franke - TPS 

28: Lawson 

19: Nielson MinNorm 

21 : Hardy Quadric 

23: Duchon - TPS 

27: Hardy Recip. Quad. 

30: Foley III 

SPLASH 

Deviations from Steep 
Table D 


.0388 

.00223 

.00447 

.0434 

.00242 

.00510 

.0317 

.00215 

.00436 

.0206 

.00176 

.00337 

.0218 

.0016? 

.00361 

.0212 

.00171 

.00337 

.0284 

.00212 

.00418 

.0210 

.00154 

.00323 

.0195 

.00101 

.00229 

.00280 

.00312 

.00031 

.0175 

.00053 

.00217 

.00736 

.00053 

.000/8 

.0143 

.00172 

.00282 

.02G3 

.OOIG 


surface, 

ICO points 

U- 



1: Franle - 3 

.143 

.0162 

.0298 

; At im a 

.115 

.0120 

.0240 

t): A Una Mod. 1 

.109 

.0113 

.0227 

13- Nielson - Fran be Q 

,0835 

.0104 

.0181 

14. Mod. Quad. Shepard 

.110 

.0113 

.0220 

16: Akin Mod. Ill 

.115 

.0119 

.0240 

,'4’ France TPS 

.150 

.0143 

.0305 

71. Lahson 

139 

.0123 

.0289 

14 Nielsen MinVom 

.115 

.0105 

.0228 

71 • Hardy Quadric 

.0716 

.00550 

.0143 

71 Pjchon - TPS 

.1 19 

.0130 

.0296 

7/’ Hard, Recip. Quad 

.0363 

.(U.r'8 

.0130 

.0 Foley Ml 

.110 

.0143 

.0249 

6SPMSU 

./AU>7 

.SI S 1 ? 


Deviations 

fren Steep test surface, 

33 points V~ 



Tabic 0.2.5 




1* 

Franlc - 3 

.113 

.0176 

.0257 

4 

At i ' a 

.0534 

. 01 23 

.0149 

in 

AU" a Mad. I 

.0320 

.0103 

.0140 

1} 

‘nelson - Fringe Q 

.0350 

.02. *w 

.0127 

u 

u oJ. Cud She; ard 

.0468 

.com 

.0126 

16 

AU" 1 Mod Ml 

.0510 

.02~23 

.0173 

« % 

Franle - TPS 

.0317 

.007* a 

.01 CD 

M. 

Lawson 

.0455 

.oi:-> 

.0255 

1 i 

fiiel'on Mnrhorn 

.0314 

.02 1 5 ’ 

.034' 1 

71. 

Hardy Q\.adi’ic 

.0138 

.024 a 3 

.00: >65 

7 3 

FucK.n - TPS 

.0233 

.02 U * 

.C05r3 

7 7 

H*rdv I cm Quad. 

.0141 

.0". 1 

. 00 ' ’ 6 

n 

F ol ry MI 

.0/13 

.010/ 

0161 



*0 l rO2. 




Deviations ft o') Strop test 

surface. 25 points 
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Generalized Cross Validation 


by Douglas Bates 
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Smoothing Surfaces 
Cross Validation 
Generalized Cross Validation 
Computational Difficulties 
Applications to Classification Problems 
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Classification Applications 


18 


_t multivariate observations 

2 populations (can be >2) 

training samples from populations 1 & 2 

Set Y-j = 1 for each t- in population 1 
= 0 for each t^ in population 2 


k 

Like noisy data generated by real valued fundtion in R 


Function is relative likelihood of sample 1 at that point. 


Fit surface using Cross Validated Multivariate Thin Plate Splines 
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to Engineering Problems 
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THE 

UNIVERSITY 
OF UTAH 


ORIGINAL PAGE IS 
OF POOR QUALITY 


DEPARTMENT Of 
MATHEMATICS 

SALT LAXE CITY IT AH 44112 
401 M1-6S51 


May 20, 1982 


Professor Larry L. Schuraker 
Department of Mathematics 
Texas A & M Uni' T ersity 
College Station, TX 77843 

Dear Larry: 

Thank you for your kind attentions at the NASA Workshop on Surface Fitting. 

You and your colleagues did a nice job of making us all feel at home. I 
regretted missing the panel discussion (and the party!). Perhaps you can 
let me know scstc of what was said. 

I have bocn thinking about what might be useful to NASA on Surface Fitting. 

The only talk that really was on NASA problems was Heydom's talk, which 
focused on statistical methods. I happen to have had a slight exposure to 
LANDSAT data from seme University of Utah geographers. Tne surface methods 
customarily used seemed to be piecewise linear or piecewise bilinear, which 
is a bit naive. Heydom's talk contained an interesting picture of fields 
of wheat, com, or "idle". This screamed for Little's arbitrary quadrilateral 
patches and/or Gregory and Cnarrot's putting triangular patches into a system 
of rectangular patches. 

A rendering issue- I was surprised that NASA thinks it can understand surface 
frcm flat pictures. I think that interactive graphics rendering is the 
absolutely bare minimum for having the illusion of understanding 3D surfaces 
A milled model is better. Some of these points were made in the DoE article 
by Barnhill and Chang which Chang referenced. This document might have some 
utility toward RASA applications - let me know if you have a copy or not. 

For arbi trarily spaced data when there is a lot of data, I think that adaptive 
methods, such as in Vittitow's PhD. thesis at Utah, are a good idea. 


In conclusion, I have the following broadly-based thoughts: 

1. The richness of possibilities for surface representation: The most im- 

portant thing about surfaces is to get one. Operation counts and all 
that pale in comparison to getting some solution to the problem at hand. 
Once a solution is found, it becomes rather routine to improve it. 


2 . 


Multidimensional problems 3D and 4D surfaces are what have resea: ch 
significance. Curves have been over-studied, even though there remain 
unanswered questions, e g. , parameterization. But just to say "3D surface" 
is insufficient. One must tailor the surface to the problem at hand. 
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3. What can NASA do with its money to achieve Surface research? 

a) Research grants to individual groups. 

b) Consultants for their labs. 

Tnese are both standard and will produce results. 

Let me suggest something new: 

c) Research grant to two groups who would collaborate on NASA- 
related Surface problems. (E.g., Barnhill's group and Schunaker' s 
group.) 

Why could this be good? Answer: the cross-fertilization would produce 

more than the sun of the parts. This would be true of both the research 
and application aspects. 

We can explore these thoughts at greater length after I hear from you. 
Cotirnunication • My office phene number is (801), 581-7916 and, if I'm not 
there, Ms. Sylvia Morris' number is (801), 581-7710. A choice of tines to 
call back would be useful. 

Best regards. 


Robert E. Barnhill 
Professor of Mathematics and 
Professor of Computer Science 

/snm 


cc and thanks: Professor Larry F. Guseman, Jr. 
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